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I.  INTRODUCTION 


During  the  past  two  years,  we  have  worked  on  several  tasks 
in  a  variety  of  areas. 

Analytical  and  computational  support  was  provided  in  the 
development  of  a  theoretical  basis  for  the  generation,  detection 
and  processing  of  microwave  magnetostatic  surface  wave  signals 
using  periodic  structures. 

A  program  was  written  to  analyze  the  diffraction  of  an  incident, 
linearly  polarized  plane  wave  at  arbitrary  incident  angle  by  one  or 
several  metal  strip  grids. 

A  program  was  written  to  evaluate  some  multiple  integrals  having 
to  do  with  the  detection  of  radar  signals  in  the  presence  of  background 
noise  and  ground  clutter. 

Programing  was  done  to  produce  a  computer- drawn  map  of  the 
earth,  or  a  certain  portion  of  the  earth. 

Programing  support  was  given  tc  ongoing  research  involving 
continuous  measurements  of  many  ionospheric  and  atmospheric 
parameters. 

Modifications  were  made  to  an  existing  program  which  had  been 
written  to  calculate  the  effects  of  different  model  ionospheres  on  VLF 
signals  reflected  from  them. 

A  computer  model  has  been  developed  to  simulate  the  performance 
of  a  ground  based  radar  so  as  to  include  the  effects  of  terrain  screening. 


bird  clutter,  ground  clutter,  Raleigh  distributed  noise,  multipath 
and  surface  roughness. 

A  program  was  written  that  computes  the  improvement  factor 
for  the  arrested  synthetic  aperture  radar. 

In  conjunction  with  the  Loran  problem,  a  program  was  written 
to  compute  the  spectrum  of  an  ensemble  of  signals. 

Finally,  analysis  and  programming  was  done  for  the  maximazation 
of  gain  of  a  wire  antenna. 
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n.  ANALYTICAL  AND  COMPUTATIONAL  TECHNIQUES  FOR  THE 

ANALYSIS  OF  PERIODIC  MAGNETOACOUSTIC  SURFACE  WAVES 

1 .  Introduction 

We  have  been  privileged  to  provide  collaborative  effort  in  the 
analytical  areas  of  numerical  and  computational  support  in  the  develop¬ 
ment  of  a  theoretical  basis  for  the  generation/detection  and  processing 
microwave  magnetostatic  surface  wave  (MSSW)  signals  using  periodic 
structures. 

This  work  describes  the  interaction  of  magnetic  surface  waves  with 
multiple  and  periodic  conductors  carrying  nonuniform  currents.  In 
addition,  a  double  ground  plane  structure  is  included  in  the  analysis, 
and  expressions  are  derived  for  radiation  resistance  and  space  harmonic 
amplitudes.  The  results  are  applicable  to  the  advanced  design  of 
magnetostatic  surface  wave  filters  and  delay  lines. 

2.  Statement  of  the  Problem 

A  transducer  in  the  form  of  a  meander  line  or  grating,  such  as  shown 
in  Fig.  1-1  is  excited  with  RF  current.  Fig.  1-2  illustrates  how  the 
transducer  is  connected  to  the  ground  plane  structure  and  to  the 
input/output  line.  The  current  sets  up  magnetic  fields  which  generate 
a  variety  of  propagating  modes  within  the  structure.  One  particular 
mode  has  a  high  potential  for  sophisticated  signal  applications 
at  microwave  frequencies.  The  mode  known  as  a  magnetostatic  surface 
wave,  MSSW,  has  the  following  characteristics:  it  is  non-reciprocal; 
propagates  perpendicular  to  a  magnetic  biasing  field  and  is  guided  by 
two  parallel  surfaces.  Most  of  its  energy  is  concentrated  near  one 
surface;  it  propagates  with  a  velocity  between  that  of  acoustic  and  light 
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GROUND  AND  SHIELD 


waves;  its  velocity  is  magnetically  tunable;  it  handles  milliwatts  o{ 
power;  and  in  practical  situations,  its  frequency  range  is  between 
2  and  20  GHZ.  The  modes  generated  by  the  structure  of  Fig.  1-1. 
are  basically  coupled  electromagnetic-ferromagnetic  modes. 

Maxwell's  equations  and  the  gyromagnetic  equation  with  electromagnetic 
boundary  conditions  must  be  simultaneously  satisfied.  The  geometry 
for  the  analysis  is  shown  in  Fig.  I- 1 .  With  reference  to  this  figure, 
the  equations  which  must  be  satisfied  by  the  fields  surrounding  the 
periodic  structure  and  bounded  by  the  two  ground  planes  are: 

-  |r-  •  *  •  « -  ° 

(D 

v  x  ft  =  ~  ^  ,  v  .  ft  =  o 


The  constitutive  relation  in  all  three  regions  are: 


ft  =  M  (ft  +  M)  =  \i  ij  .ft 

o  o  y 

5  =  .  t 


(2) 


In  regions  (I)  and  (III),  M  =  0,  and  the  gyromagnetic  equation  for 
region  (H) is: 


at 


-  _  Y  M  x  ft 


(3) 


This  gyromagnetic  equation  is  to  be  linearized  to  first  order  in  small 
signal  RF  field  variables. 
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3 .  T r ansverse  Electric  (TE)  Modes 

Maxwell's  equations  can  be  separated  into  two  sets  of  equations: 

one  set  for  E  ,  E  ,  and  H  (the  transverse  magnetic  or  TM  mode), 
x  y  z 

and  one  set  for  E  ,  H  ,  and  H  (the  TE  mode).  TE  mode  solutions 
z  x  y 

are  characterized  by  the  following  conditions:  B/dz  -»  0  for  all  fields  and 
the  wave  solutions  have  the  form 

fk  (y)  exp  [i  (uJt  -  kx)] 

where  fk  (y)  =  exp  (ky)  +  BkR  exp  (-ky) 

for  every  region  R  (Fig.  1-1).  The  magnetic  biasing  field  H  is  Z 
directed.  The  nonzero  field  components  satisfy  the  following  reduced 
equations. 


dHy 

dH  _ 

-  x  =  uufiE. 

Z 

(4) 

B* 

ay 

*Ez 

=  -  uuB 

ay 

*Ez 

=  io)B 

y 

(5) 

a* 

*BX 

4  8By  =  0 

ax 

oy 

The  V  •  D  equation  is  automatically  satisfied  since  E  =  E  =0  and 

x  y 

d/d z  -*  0.  The  magnetization  vector  has  the  form 


M.  =  M  X  +  MY  +  MZ 
x  y  o 
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and  the  linearized  magnetization  equation  reduces  to: 


i  ujM  =  vM  H  -  VH  M 
x  o  y  1  o  y 

i  uiM  =  yM  H  4  yH  Jvd 
y  ox  'ox 

Using  ^  4  M)  and  the  preceding  equations. 


(7) 


N  =  i 

VlZ\ 

i H  \ 

IbJ  °( 

^'lu12  u22/ 

(8) 


The  solutions  for  B  and  H  in  each  of  three  regions  are  given  as 

y 

integrals  of  functions  f^y)  =  A^exp  (ky)  4  BkRexp(-ky),  where  R 
is  the  number  of  the  region. 

The  electric  field  in  each  of  the  three  regions  is  obtained  by 

multiplying  the  integrand  of  B  by  (-tu/k).  H  and  B  are  obtained 

y  y  x 

from  B  and  H  from  the  permeability  matrix  (8). 

y  * 

All  of  the  field  equations  and  the  final  power  expression  are  valid 
for  anisotropic  media.  For  the  isotropic  case 


nH 

Wn  '  M22  "  1  '  tFTi?" 

H 

=  “z  Tz~  ’  =  H/(4nM0)  (9) 

H 

Y  =  2.  8  MHZ/oersted 

a  _  UL  j 

y  .  4 rtM  4tcM  =  1750  oersted 

o  o 
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4.  Boundary  Conditions 

In  order  to  solve  the  equations  of  the  preceding  section  and 
evaluate  the  constants  A^>  B  boundary  conditions 

apply.  With  reference  to  Fig.  1-1, 


B 

y 


0  at  y  =  t  j  and  y  = 
B  11  at  y  =  -d 

y 


(d  +  L) 


H 


II 


at  y  =  -d 


(10) 


The  boundary  condition  due  to  current  in  the  periodic  structure  is 


H  m  -  H  11  =  J(x)  at  y  =  0 

X  X 


(11) 


All  coefficients  are  given  in  terms  of  a  Fourier  transform  of  an 
arbitrary  surface  current  distribution  J(x). 

This  is  an  important  result  because  it  provides  a  procedure  for 
sophisticated  transducer  designs.  For  example,  non-periodic  (or 
periodic)  conducting  strips  may  be  driven  with  independent  current 
sources;  then  an  integration  (or  Fourier  transform)  over  the  source 
distribution  provides  a  quick  determination  of  all  important  trans¬ 
ducer  characteristics. 


5 .  General  Fields  and  the  Characteristic  Equation 

To  obtain  the  total  fields  in  the  three  regions,  which  are  integrals, 
use  of  residue  theory  is  required. 
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The  characteristic  relationship  between  frequency  and  wave  number 
for  the  unelectroded  structure,  consisting  of  Y1G  slab  and  both 
ground  planes,  is  given  by 


-2fB2  |k|d 


FT(k,  tu)  =  P3(a2e  -  OfjT)  coth  (P3  Jk  |tj) 

-2p2/k/d^  T[tfip3COth(P3|k|t1)  +1] 


-2e_  |k|d 


+  T)  =  0, 


where 


a2P3coth(p3  |k  (t^. 


2  =  - 
^1 1^2ivl12  S 

T  =  °2  +  UPitanl1  (Bj  |k|t) 

Q-J  -  1/3  jtanh  (gj  I k  |  "-) 


S  =  ±1. 


6.  Magnetostatic  Surface  Solutions 

In  the  magnetostatic  limit,  MSSW  wavelengths  are  much  smaller 

than  electromagnetic  wavelengths  in  regions  of  dielectric  constant  e. 

That  is,  X  si  .  This  is  equivalent  to 
s  e 


k  »  Jt  uj/c  -  k 
s  o  e 

9  1  8 

UJ  =  2ttx  3  X  10  sec  ,  cq  =  3  x  10  m/sec  , 

o  -  1  o 

£  =  10.  It  was  found  that  k  «  2  cm  or  X  TT  cm. 

C  £ 


This  means  MSSW  wavelengths  up  to  a  few  thousand  microns  satisfy 
inequality  (13)  very  well.  Practical  MSSW  wavelengths  are  on  the 
order  of  hundreds  of  microns  and  so  (13)  is  well  satisfied  in  practice. 

For  a  given  frequency  u),  k  is  found  from  the  characteristic  equation 


F  =  F 
TM  T 


k»k 


o  =  (o^e-2pMd  -  fflTM)  coth  (  |k  Jtj)  -  (e'Zplk  id  +  TM)  , 


Where  3  = 


(14) 


T 


M 


+  tanh  (  |  k  1 1) 
-  tanh  (  |k  |-l) 


7.  MSSW  Fower 

From  the  V  x  2  =  -2 


equation, 


E  (S)  =  -S(u>/k)  B  (S) 
z  y 


(15) 


and  from  the  permeability  matrix,  (8) 


„  (s)  .  /B  (s)+  in  ?H  (S)' 

Hy  =  lf»2Z  _  12  * 


(16) 


The  power  density  carried  by  a  MSSW  in  the  x  direction  is 


P(s)  =  1/2  E  (S)H  (S) 

z  y 


(17) 
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while  total  power. 


for  width  t  , 


is 


E  (S)  .  H  (6)dy  (18) 

z  y  7 

-(■*+  d) 

After  substitution  of  B  H  ^F_w,  E  and  H  ^  into  P^,  a 

y  x  TM  z  y 

final  expression  for  P  is 


,(s)  =  ^oS  |c(k,  o))|  2|  A(ks)j, 
2k  2 


(19) 


Where 


A(ks)  =  (TM+  1) 


Sinh  (2k  l)  k  L 
s  -  s 


,-pk  d) 


-  «2e  r  s 


Cosh  (k  -t.) 
Sinh  (2kgtl)  -  kgtj 


Sinh2  (k  t.) 

si 


-2Sk  d 


a2  (e 

T 


-  >>  -  2ksdTM^ll 


,  „  Bk  d 

4  (orjT^e  s 


_  2Rk  d 
,  O'.  rp  2  6 

+  -1  tm  <e  -D- 


ln  evaluating  power,  (1  9),  k^  =  |k|  satisfying  =  0  from  (14),  and 


G(k, 


-Bk  d 

«,)  -  7(5  •  V  6 _ ^ 

FTM  (ks’  tu) 


(20) 
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8.  Radiation  Resistance 


In  this  section  are  given  expressions  for  radiation  resistance  of 
a  meander  and  grating  line  transducer  assuming  uniform  current 
distribution  in  each  conductor. 

A  transducer  is  made  up  of  N  conducting  strips  each  of  width  ”a”  carrying 
current  I  . 


When  T)  =  +1  all  conductors  are  connected  in  parallel  to  form  a  grating, 
and  when  7)  =  -  1  they  are  connected  in  series  to  form  a  meander  line. 

After  taking  the  Fourier  transform  of  this  current  distribution  and  substitution 

/g) 

into  (19)  and  (20), the  radiation  resistance,  R  ,  can  be  obtained  as 

m 

follows: 

p,s)  -  >/*  •  M 2  ’  <21) 


Where  I  =  \JZ  [(1  -  T))  +  (1  +  7])  N] 


[  (1-  T))  +  (1  +  T))  N2] 


f  Sin  ( 

L  W 


a  k  JZ) 

W 


I  -  T)NeikPN 
1  -  7)eikP 


(s)  ^oV 

"  ,  .2  lr  (1) 


-2gk  d 
s 


ks  |FTM  <V  H 


Then  for  a  grating,  T(  =  +  1, 

.  .  R  (S)  Sin  (ak  /2) 
p  _  _£>  s 

Km  '  (akg/2) 


Sin  (k  PN/2) 
_ s _ 

Sin  (kgp/2 


For  a  meander  line,  7)  =  -1,  with  N  even. 


R  (,)  .  i.)  Si"  (a  V2) 

m  o  (a  k  JZ) 

s 


Sin  (kgp  N/2) 

Cos  (k  P/2) 
s 
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Fig.  1-3  gives  curves  of  total  radiation  resistance  per  unit  width, 
(Rm  =  Rm  f°r  the  single  strip  case,  N  =  1. 

p  (+)  .  _  (-) 

m  and  Rm  are  obtauled  using  equation  (22).  For  all  curves, 
a  =  178  |im,  H  =  650  oe  and  d  =  6.  25  pm. 


Computational  Techniques  and  Numerical  Analysis 


Program  EXC  has  been  developed  to  effect  numerical  analysis 
m  order  to  test  the  general  design  theory  for  magnetostatic  surface 
wave  periodic  transducer  structures;  analyze  the  behavior  of  the 
functions  P  and  R;  obtain  an  accurate  correspondence  between  theory 
and  experiment;  and  compare  the  results  thus  obtained  with  those  of 
other  investigators.  This  testing  and  analysis  consists  of  performing 
the  following  functions: 

Plotting  Dispersion  Diagram  for  0  =  f(k)  for  both  F  (k,  w)  =  0  (S  =  ± 
and  FTM(k’  “•')  =  0  (S  =  ±  1). 

Calculation  of  the  power  carried  by  MSSW  P^  (S  =  1)  and  P  (S  =  -1). 

Analysis  of  the  behavior  of  the  power  and  Radiation  Resistance 
functions  in  the  neighborhood  of  a  given  frequency  a«  and  plotting 
the  graphs  of  these  functions. 

Evaluation  of  the  power  generated  as  a  function  of  width  to  half 
period  ratio  (a/p). 

To  carry  out  the  functions  mentioned  above.  Program  EXC  calculates  the 
roots  of  transcendental  equations  F Tjvl(k,  (u)  =  0  and  F  (k,  a-)  =  0  for 
every  given  frequency  uj  relating  wave  number  k. 


FREQUENCY 


The  main  computational  problem  was  to  select  a  suitable  method  for 
solving  transcendental  equations. 

Existing  methods  for  the  numerical  solution  of  transcendental 
equations  generally  apply  to  continuous  functions  with  continuous 
derivatives.  Furthermore,  nearly  all  of  them  require  some  prior 
knowledge  of  the  neighborhoods  of  the  roots.  Because  the  functions 
associated  with  the  above  analysis  are  not  in  general  well-behaved, 
and  because  it  is  sometimes  impossible  to  predict  all  the  root 

locations,  it  has  been  necessary  to  make  appropriate  modifications 

6 

to  the  Mueller  iteration  scheme.  This  was  done  by  insertion  of  the 
condition  at  infinity  for  transcendental  functions  within  the  given 
range  of  argument. 

On  the  other  hand,  subroutine  RTMI  requires  knowledge  of  the 

X . .  and  X  .  left  and  right  bounds  of  the  range  of  X,  where  unknown 
*01  ri 

roots  exists.  To  overcome  this  restriction,  the  whole  range  of 

frequency  (or  wave  number)  is  divided  into  n  equal  segments,  for 

each  of  which,  DISCR  =  X..  •  X  .  is  calculated.  That  pair  of  X  and 

li  ri  r  li 

X^.,  which  gives  a  negative  value  of  DISCR,  are  left  and  right  bounds 
for  the  root.  For  oscillating  functions,  n  should  be  large 
(from  100  to  300). 

A.  Plotting  Dispersion  Diagram 

The  Dispersion  Diagram  was  plotted.  For  each  value  of 
frequency  uu,  the  equations  FT(k,  ou)  =  0,  FTM  (k,  uj)  =  0  were  solved 
relating  wave  number  k  for  both  cases  (S  =  ±  1). 
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B.  Calculation  of  the  Power  Carried  by  MSSW 

The  calculation  of  the  power  carried  by  MSSW  (P  and  P  )  was 
realized  by  means  of  subroutine  POW,  varying  frequency  u>  for  given  S. 

The  value  k,  which  occurs  in  the  function  POW,  is  calculated  using 
subroutines  RTMI  and  FTM. 

C.  Plotting  the  Graph  of  the  Radiation  Resistance 

Plots  of  the  Radiation  Resistance  as  a  function  of  frequency  were 
performed  on  the  CDC6600  on-line  plotter.  The  frequency  range  is  given, 
the  number  of  arguments  varies  depending  on  the  goal  of  the  analysis; 
k-values  are  the  same  as  for  power  calculation.  The  range  of  frequency 
is  not  divided  into  equal  segments  because  in  the  neighborhood  of  the  peak, 
(eg  =  3.56  GHZ)  the  function  is  rapidly  increasing  and  more  points  are 
needed. 

D.  Evaluation  of  the  Power  Generated  as  a  Function  of  the  Parameters 
The  power  (or  Radiation  Resistance)  is  very  sensitive  to  variation  in 

the  parameters.  Variation  of  the  parameters  a,  L,  t^  must  be  in  the 

range  over  which  the  function  exists  (Fig.  1-3).  For  example,  with  the 

parameters  given  in  the  Figure,  the  junction  (Radiation  Resistance) 

9 

does  not  exist  below  about  22  x  10  radians/sec.  or  3.5  GH  . 

z 
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IU-  analysis  of  propagation  through  multiple  layers 

OF  METALLIC  GRIDS  (PROGRAM  ALT) 

Program  ALT  has  been  written  to  analyze  the  diffraction  of  an 
incident,  lineary  polarized  plane  wave  at  arbitrary  incident  angle  by 
one  or  several  metal  strip  grids. 

Each  strip  has  zero  thickness, and  width  a  (for  the  g-th  grid). 

The  grids  each  have  the  same  spacing  d^  and  are  located  at  positions 
Zg.  Each  strip  is  infinitely  long,  and  there  are  an  infinite  number  of 
strips  in  each  grid. 

The  incident  wave  has  an  arbitrary  angle  of  incidence  along  the 
radial  vector  in  conventional  polar  co-ordinates  (r,  0,  cp),  with  0 
measured  from  the  z-axis  and  Cp  measured  from  the  x-axis  of  a 
rectangular  co-ordinate  system.  The  incident  wave  has  its  E-field 
polarized  along  the  plane  formed  by  z  and  the  angle  of  incidence. 

Both  components  of  scattered  field  are  included.  The  parameters 
u  and  v  are  the  direction  cosines.  These  are  inputed  to  the  program 
as  u  and  v.  is  the  number  of  layers.  The  analysis  assumes  an 

indire ctional  current  along  the  Y  direction  in  each  strip  and  solves 
for  the  complete  diffraction  field  including  primary  and  cros s- polarize 
field  and  the  radiation  into  grating  lobes. 


1  •  Computational  Techniques  arid  Numerical  Analysis 

The  purpose  of  program  ALT  is  to  compute  and  print  out  the  five 
components  of  the  total  field  along  with  the  values  20  logio  of  their 
magnitude . 


These  components  allow  the  analysis  of  the  propagation  effects 
through  a  filter  consisting  of  thin  metal  strips  spaced  d  apart,  with 
their  axis  oriented  along  the  same  direction. 


The  input  parameters: 

°M'  {am},  {zm}'  dx’  U'  V’  (m  =  1,  2,  ... 


The  output  (components  of  the  total  field): 

tjp),  r(P),  T(P),  c  +  (P).  c  -  (P) 

TsdB(P)'  rdB(P)'  TdB(P)'  C+  dB(P)'  C-dB(P) 

Y,  CONTR  (control  sum  for  the  checking  of  calculation). 

These  components  are  to  be  printed  for  each  value  of  P,  subject  to 
the  condition  described  below. 


The  complete  solution  breaks  up  into  four  parts: 


1.  Finding  the  elements  of  the  vector 


where  e 


m 


m 


m 
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2.  Finding  the  elements  (complex)  of  the  matrix  [c] 

•where 


3.  Solve  the  Matrix  Equation: 

W  \c\  =  {j} 


4.  For  each  p,  subject  to  the  condition: 


29 


compute: 


t.(p»  =  F+  (p)  k0  kP  Y=ffm 

r  (P)  =  f_  (P)  ko  kp 

T  (P)  =6  (P)  +  Tg  <P) 

C  +  (P)  =  -  F+  (P)  Pp  ky 
C-  (P)  =  -F_  (P)  p^ky 


along  with  20  times  the  logarithm  of  their  magnitude, 
where 


6(P)  =  1 

F+(P)  = 

F_(P)  = 


for  p  =  0 

_ 1 _ 

2  2 

Z(k  -ft  )  k  d 
o  Kp  ox 


-jk  z 


E  jq  — - 9  sin^paqj 


-2(k  2  -B  2)  k  d  9 
o  'ox 


jk  z 

E  jq  -  -  ■ ^ — 9  sin^p3 
bp 


The  control  sum  of  the  solution  is: 


P 


=  Z>[  J  T(p)  I2  +  I  r(P)  I  c+(P)  |2+  I  c_(P)  |2)]  ~b 

o 


and  should  be  equal  to  unity. 
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Because  the  elements  of  the  matrix  [c]  are  an  infinite  series  of  complex 
terms,  the  left  and  right  bound  of  the  sum-PM  are  chosen  from  1  to  20, 
depending  on  the  value  of  input  parameters.  (This  series  converges  very 

rapidly).  In  the  process  of  the  solution  of  the  matrix  equation,  the  sub- 

2 

routine  MINV  ,  modified  for  complex  elements,  was  used. 


2.  The  Input  and  Output  of  Program  ALT 

The  input  to  program  ALT  is  contained  in  namelist  INP.  The 
parameters  of  namelist  INP  are: 


QM 

-  size  of  the  matrix  [c] 

A  (I) 

-  value  of  the  I-th  element  of  the  vector  ia  ) 

l  mj 

f 

Z(I) 

-  value  of  the  I-th  element  of  the  vector  / z  1 

\  ml 

1 

DX 

-  value  of  the  d 

X 

V 

-  value  of  the  v 

U 

-  value  of  the  u 

PM 

-  value  of  the  left  and  right  bounds  of  the  sum  (for  the 

calculation  elements  of  the  matrix  [ c] ) . 

31 


The  output  of  program  ALT  is  printed.  The  printed  output  is: 


1)  The  namelist  INP  (described  in  input). 

2)  The  namelist  OUT.  The  parameters  of  namelist  OUT  are: 

P  -  value  of  p  subject  to  the  conditions  described  in 

FUNCTIONAL  DESCRIPTION. 

T  -  T(P)  -  value 


TS 

-  T  (P)  -  value 
s 

GAM 

-  r(P)  -  value 

CPLUS 

-  C+(P)  -  value 

CMIN 

-  C  (P)  -  value 

TSDB 

-  T  value 

sdB 

TDB 

-  TJn  -  value 
dB 

GAMDB 

-  r  -  value 

dB 

CPLUSD 

'  C+dB  '  ValUC 

CMINDB 

-  C  -  value 

-dB 

Y 

-  Y  value 

CONTR 

-  P-  value 

If  any  of  the  values  don't  exist,  the  corresponding  message  is 
printed. 

The  value  of  CONTR  (on  output)  should  be  equal  to  unity.  If 
CONTR  is  not  equal  to  unity,  it  means  that  input  data  (parameters) 
don't  satisfy  the  restrictions  1-4. 
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IV.  DETECTION  OF  SIGNALS  IN  NOISE  AND  CLUTTER 
1 .  Introduction 

The  request  was  made,  to  numerically  evaluate  some  multiple 

integrals  having  to  do  with  the  detection  of  radar  signals  in  the 

presence  of  background  noise  and  ground  clutter.  In  particular,  the 

integrals  P  and  P  ,  given  below,  are  probabilities  that  a  radar 
dl  dZ 

detection  system  identifies  an  existing  target,  such  as  a  missile  or 
airplane,  while  set  at  a  particular  threshold  level  and  receiving  N 
reflected  pulses  per  sweep.  The  integral  is  the  probability  that 
a  "false  alarm"  occurs:  that  a  "target"  be  seen  when  none  is  actually 
present. 

The  difference  between  P  and  P  is  that  in  the  latter  case, 

dl  dZ 

any  fluctuations  in  the  target  are  assumed  to  occur  slowly  in  compari¬ 
son  to  the  repetition  rate  of  the  pulses  within  each  sweep,  whereas  in 
the  former,  the  fluctuations  may  be  rapid.  Thus,  for  the  special  case 

N  =  1,  the  triple  integral  P,,  must  reduce  to  the  double  integral  P,_. 

dl  dZ 

This  was  checked  analytic  ally.  Moreover,  for  the  special  case  in 
which  there  is  no  signal  (signal  to  noise  ratio  zero),  P^  reduces  to  P^ 
regardless  of  N.  The  assumptions  are  made  that  the  signal  and  the 
noise  are  Rayleigh- distributed,  and  that  the  clutter  has  a  long-normal 
distribution. 

The  integrals,  transformed  to  the  most  convenient  form  for  compu¬ 
tation,  are 


Pdl  =  1  * 


=”J  e"ydyjd§ 


-(1  +  Na)§ 


zjz  no  o 


GN(y,Na5)  J  ds -e  S yE^^sT)  * 
o 

r  1  2 

exp  J  -  — -In  (2q  (3  s)  V 
8q 


o 

J  e'ydy  f ds 7 e_s  cN,y- s)  exp  {-  ?i2 

tt  0+  r\  n  u  ' 


Pd2  =  1 


zjzvc  o 
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g 


is  identical  to  P  when  the  parameter  a  is  taken  to  be  zero.  The 
GN(y,  s)  are 


N 


r-N-1 

In_1(  2,/ys) 


I^j(x)  is  the  Modified  Bessel  Function  of  order  N  and  argument  x, 

N  =  number  of  pulses  per  sweep, 
a  =  signal  to  noise  ratio, 

P  =  noise  to  clutter  ratio, 

a  =  variance  of  the  log-normal  clutter  distribution, 

=  threshold  level  setting 

The  objective  (for  each  desired  N,  g  and  a)  was  to  calculate  Pf  as 

a  function  of  y^,  thereby  establishing  values  of  y  for  which  P^.  was 

satisfactorily  small;  and  then,  using  these  y  ,  to  find  P,.  and  P  as 

o  al  a2 

a  function  of  a,  proceeding  from  small  a  to  values  large  enough  that 
the  probability  of  detection  reaches  at  least  .  99. 

2.  Results 

The  parameters  used  are  listed  below: 

N  =1,3,  6,  10 


3  =  J. T,  1,  ,/n),  10,  100 

0  =  0. 7,  1. 1,  1.6 
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The  values  of  y  required  to  give  sufficiently  small  have  the 

following  dependence  on  the  parameters:  the  larger  p  is,  the  smaller 

y  ;  the  larger  N  is,  the  larger  y^;  and  the  larger  o  is,  the  larger  y^. 

The  dependence  on  N  and  p  is  strong,  but  the  dependence  on  o  is  even 

more  so.  In  fact,  values  of  y  for  o  =  1 .  1  and  1 . 6  were  so  large  that 

o 

the  computer  time  required  to  do  the  integrals  was  considerable,  and 
these  cases  were  dropped  after  a  cursory  investigation. 

Table  1  gives  the  values  determined  for  y^  for  different  cases. 

By  making  the  assumption  that  y^  >>  1,  one  can  show  that 


where  erfc(x)  is  the  complementary  error  function.  This  expression 
is  quite  accurate  for  p  £  J\0,  and  serves  as  a  check  on  the  numerical 
results.  As  p  increases  above  5,  the  approximate  expression  becomes 
progressively  less  accurate. 

P,_  was  then  calculated  as  a  function  of  a .  The  function  rises 
ad 

sharply  (from  f^.)  as  a  increases  from  zero,  then  flattens  out  and 
slowly  approaches  1  as  a  becomes  large.  This  is  illustrated  in 
Figure  1  for  a  typical  case.  The  larger  p  is,  the  larger  P  is  for  a 
given  a •  (Equivalently,  for  larger  p,  smaller  a  are  required  to  get  a 
fixed  detection  probability.  )  In  the  flat  part  of  the  curve,  larger  N 
give  larger  values  of  P  . 

It  was  possible  to  derive  a  limiting  expression  for  P  which  is 

ad 

valid  in  the  limit  of  large  a : 
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This  provides  a  quick  and  accurate  answer  for  a>y  and  also  serves 
as  a  check  on  the  numerical  results  for  such  cases. 

Tables  2  and  3  give  the  values  of  a  for  which  P  is  equal  to  .9, 

4  a  ^  ^ 

.  99,  and  .999,  when  =  10  and  P^  =  10  ,  respectively. 

P.,  was  also  calculated  for  different  a •  For  N  =  1,  the  results 
di 

should  be,  and  are,  identical  to  F  .  For  N  >  1,  the  function  behaves 

a  2 

much  the  way  P  does  when  R  is  varied.  However,  a  variation  of  N 
a<! 

has  very  little  effect:  the  results  are  similar  for  all  N. 

A  limiting  expression  was  derived  for  P  for  large  a: 


It  is  possible  to  do  the  §-  and  y-  integrals  in  P  analytically.  The 
result  is  an  infinite  series,  with  Gauss  hypergeometric  functions^  in 
the  integrand  of  the  s-integra!  in  each  term.  The  convergence  proper¬ 
ties  of  the  series  are  unknown,  and  this  approach  was  not  pursued 
vigorously. 


3.  Programs 

Pj  and  P^^  are  calculated  in  program  PPWN08,  and  P^  is  calcu¬ 
lated  in  PPW12.  The  input  data  for  each  program  are  read  from  a 
single  card  for  each  run.  These  are  described  in  detail  in  the  program 


4.  PPWN08 


The  repeated  computation  of  the  s-  integral  is  the  most  time- 

consuming  portion  of  the  job.  The  front  edge  of  the  s-  integral 

(0  £  s  <  1,  where  the  integrand  varies  most  rapidly)  is  done  with 

2 

an  adaptive  Simpson's  rule  routine.  This  uses  an  interval -halving 

procedure  which  minimizes  the  number  of  integrand  evaluations. 

The  rest  of  the  s-  integral,  1  £  s  s  s  ,  is  done  by  m-point  Gauss- 

2  m  s _x 

Legendre  quadrature,  where  m  is  read  in.  The  region  1  to  s 

max 

is  divided  into  pieces  (up  to  8  of  them)  whose  width  depends  on  y. 

The  integration  proceeds  piece  to  piece  in  the  direction  of  increasing 
s  until  one  piece  gives  a  negligible  contribution  to  the  total.  If  the 


8th  interval  still  contributes,  a  9th  interval  (s  to  ®)  is  processed 

max  e 

by  making  the  transformation  x  =  1/s  and  integrating  0  s  x  fi  1/s 

max 


Otherwise,  s  is  a  substitute  for  <=. 
max 


When  y  is  small,  the  s-  integrand  drops  off  slowly  as  s  increases 
beyond  y.  For  large  y,  it  is  sharply  peaked  at  s  ~  y.  In  many  such 
cases  the  region  0  S  s  S  1  contributes  negligibly,  so  the  worst  part  of 
the  computation  can  be  avoided. 


The  y-  integration  is  done  by  m' -point  Gaussian  quadrature,  m' 
being  read  in.  For  N  =  1,  the  integrand  decreases  monotonically  from 

y  =  0;  for  N>1,  it  peaks  at  small  y.  In  the  latter  case,  two  intervals 
are  integrated  separately,  for  improved  accuracy. 

The  main  program  reads-  in  the  input  data,  prints  it,  and  calculates 

the  limiting  case  (for  large  a,  to  compare  with  P,_).  Subroutine  CHOICE 

dZ  4 

chooses  the  desired  m'-point  Gaussian  quadrature  subroutine  DQGm. 
Function  FN  is  the  y-  integrand.  It  calls  CHOIC,  which  selects  DQGm, 
which  integrates  function  SINT,  the  s-integrand. 
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5.  PPW12 


The  triple  integral  P  was  transformed  so  that  the  s-  integral 

dl 

would  look  exactly  like  that  of  ,  except  that  I^^^gs)  is  substituted 
for  G^(y,  §).  The  same  comments  about  the  shape  of  the  s-  integrand 
apply.  The  region  0  S  s  S®is  divided  into  pieces  exactly  as  before, 
but  an  adaptive  4-point  Gaussian  integration  subroutine  AQG4  was 
written  to  replace  the  adaptive  Simpson's  (ASQ)  and  standard  Gaussian 
routines.  The  superior  accuracy  of  the  Gauss- Legendre  formula  off¬ 
sets  the  inherent  efficiency  (not  recomputing  any  integrands)  of  ASG 
and  provides  a  somewhat  (up  to  20%)  faster  evaluation  of  the  s-  integral, 
for  a  given  required  accuracy. 


The  integrand,  typically,  is  very  sharply  peaked  at  §  =  y/No. 
AQG4  is  used  for  this  integral  also,  as  it  is  for  the  y-  integral,  which 
is  much  like  the  y-  integral  in  P  . 

The  s-  integral,  being  the  most  deeply  imbedded,  is  required  so 
frequently  that  it  is  efficient  to  compute  a  table  consisting  of  §  and  the 
s-  integral,  and  interpolate  whenever  a  value  of  the  latter  is  needed. 
This  saves  a  great  deal  of  time. 


The  main  program  reads  the  input  data,  prints  it,  and  calculates 
the  interpolation  table  by  calling  SINTGL.  This  subroutine  performs 
the  s-  integral  by  calling  AQG4  to  integrate  function  SINT.  The  main 
program  then  calculates  a  divided  difference  table  from  the  interpo¬ 
lation  table,  to  estimate  the  worst  errors  likely  to  be  made  by  the 
Aitken- Lag  range  interpolation^  procedure  using  INT  table  entries. 

INT  is  automatically  adjusted  upwards  (to  a  limit  of  9)  if  greater 
accuracy  is  required.  It  then  does  the  remaining  integrations  by  call¬ 
ing  AQG4W,  which  integrates  function  YINT,  the  y-  integrand.  YINT 
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sets  the  limits  for  the  §-  integration  and  calls  AQG4P  to  perform  it. 
Function  CSINT  calculates  the  integrand.  Part  of  the  latter  pro¬ 
cedure  consists  of  calling  subroutine  INTERP,  which  orders  the  table 
entries,  and  ALI,  which  performs  the  interpolation  itself. 

Creating  the  table  requires  the  most  time,  so  it  is  advantageous 
to  use  it  for  more  than  one  case  at  a  time.  If  the  product  p-  a  Temains 
unchanged  from  one  case  to  the  next,  the  table  is  not  recomputed. 

The  required  range  for  the  table  (Os^s  CSIMAX)  is  slightly  greater 
for  smaller  N  for  a  given  p»  a >  so  a  case  with  N  =  1  should  be  followed 
rather  than  preceded  by  cases  with  N>1.  If  this  is  not  done  and  a 
point  outside  the  table  range  is  required,  CSINT  calls  SINTGL  to 
calculate  it  rather  than  extrapolating  from  the  table. 


6.  IN  and  INASYM 

These  subroutines  were  written  to  calculate  the  modified  Bessel 

functions  I  (x)  to  any  desired  accuracy.  IN  uses  the  infinite  series 
7 

representation  which  converges  rapidly  for  small  x,  while  INASYM 

7 

uses  an  asymptotic  series  and  is  useful  for  large  x.  Tests  show  that 
for  a  desired  accuracy  of  10  ,  IN  is  more  efficient  for  x<k  +  n  and 

INASYM  should  be  used  otherwise.  For  n  =  0  or  1,  the  cutoff  point  is 
x<k  +  2. 
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f 


N 

L 

0  = 

1 

Jio 

10 

100 

JlQ 

^10 

1 

4 

287 

94.96 

30 

13.  67 

9.35 

565 

23,265 

3 

4 

835 

279.  41 

90 

33.46 

14.  23 

1680 

6 

4 

1735 

556. 11 

180 

63.  68 

20.  19 

3400 

lu 

4 

2900 

925.  07 

300 

104. 15 

27.41 

5640 

, 

1 

t 

5 

630 

200 

66 

23 

1875 

3 

5 

1880 

5  98 

192 

62 

5650 

6 

5 

3720 

1180 

400 

129 

11294 

10 

5 

6220 

1970 

630 

208 

18781 

1 

6 

1240 

390 

125 

43.44 

14.11 

5600 

3 

6 

3700 

1170 

375 

122.  96 

20.76 

16500 

6 

6 

7400 

2350 

730 

242.42 

33.  67 

32850 

10 

6 

12400 

3880 

1240 

402.  70 

53. 00 

55860 

1 

7 

2300 

700 

235 

77 

3 

7 

6800 

2170 

700 

223 

6 

7 

4350 

1400 

440 

10 

7 

7250 

2250 

765 

Table  1.  Values  of  y^,  determined  from  program  PPWN08,  for 
which  Pj.  =  10  .  o  =  0.7  except  for  the  sixth  column,  where 


a  =  1.1,  and  the  last,  where  g  =  1.6.  Values  without  decimals 
are  accurate  to  5%;  others  are  accurate  to  the  last  decimal  given. 
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N 

Pd2 

P  =/T 

1 

yro 

10 

100 

1 

.9 

2723 

872 

284 

129 

87.  8 

1 

•  99 

2.  86  x  10 

9150 

2985 

1360 

930 

1 

.999 

2.  87  x  10 

9. 20  x  10 

3.  00  x  10 

1. 37  x  104 

9348 

3 

•  9 

750 

250 

80 

29.  5 

11.8 

3 

-99 

1900 

630 

205 

76.  0 

31.7 

3 

.  999 

4350 

1430 

470 

175 

74.0 

6 

.  9 

547 

175 

56.  0 

19.  2 

5. 47 

6 

.99 

965 

227 

99.6 

34 

10.  3 

6 

.  999 

1560 

320 

169.  6 

56  . 

17.  0 

10 

•  9 

466  . 

151 

47 

15.7 

3.4 

10 

•  99 

693 

228 

71.6 

24 

5.6 

10 

.  999 

963 

320 

100 

34 

8.  2 

Table  2.  Values  of  a,  determined  from  program  PPWN08, 
for  which  P =  .  9,  .  99,  and  .  999-  q  =  0.7  and  P  =  10  4. 
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N 

P^2 

,  p  = 

l 

yr0 

10 

100 

1 

•  9 

X 

00 

rH 

105 

3.  70  x  10^ 

1. 19  x  10^ 

411 

133 

1 

.  99 

1 .  23  x 

106 

3.  88  x  104 

1. 24  x  104 

4321 

1403 

1 

.  999 

1. 24  x 

iofe 

3.  90  x  105 

1 . 24  x  105 

4. 34  x  10 

1.41  x  10 

3 

.  9 

3380 

1055 

348 

113 

17.  7 

3 

.  99 

8500 

2680 

855 

280 

46 

3 

.  999 

1. 93  x 

io4 

6100 

1950 

640 

107 

6 

•  9 

2350 

640 

234 

75.  5 

9-  7 

6 

•  99 

4150 

1300 

412 

135 

17.  8 

6 

.  999 

6650 

2120 

660 

240 

29.  5 

10 

.  9 

1980 

620 

197 

64 

7.  6 

10 

.  99 

3000 

937 

298 

96 

11.  8 

10 

.  999 

4350 

1300 

415 

137 

16.  9 

Table  3.  Values  of  a,  determined  from  program  PPWN08, 
for  which  P  =  .  9,  •  99,  and  .  999-  a  =  0.  7  and  P  =  10'6. 
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V.  MAPPING  PROGRAM 


1 .  Introduction 

The  request  was  made,  to  modify  and  supplement  certain  subroutines 
provided  by  the  initiator,  write  some  new  ones,  and  then  combine  them  so 
that  all  the  results  would  be  available  simultaneously.  The  objective  was 
to  produce  a  computer- drawn  map  of  the  earth,  or  certain  portions  of  the 
earth,  with  any  desired  great  circle,  or  any  other  arc  or  line,  superposed 
on  it,  and  with  certain  labels  in  the  appropriate  places. 

The  complete  program  consists  of  a  main  program,  MAP,  two  sets  of 
subroutines,  and  a  block  data  program.  The  subroutines  DEC,  CONV,  and 
POINTS,  and  block  data  DIPEQ,  are  used  for  producing  sets  of  points 
(latitude,  longitude)  to  be  plotted  on  the  map.  The  subroutines  SUPMAP, 
MAPLOT,  SUPCON,  and  FRSTPT  aTe  used  to  read  these  data,  transform 
the  (lat.  Ion)  coordinates  to  a  form  (u,  v)  suitable  for  the  desired  projection, 
and  do  the  actual  plotting.  (For  pen  plotting,  subroutine  CODE  also  prints 
the  map  parameters  next  to  the  map. )  The  standard  CRT  or  PEN  routines 
are  required,  as  well. 


2.  Description  of  the  Program 
2.  1  General 

In  their  original  form,  SUPMAP,  MAPLOT,  SUPCON  and  FRSTPT  are 
capable  of  producing  the  map  itself,  which  shows  the  land  contours  with  a 
resolution  of  about  .  1°.  The  use  of  these  subroutines,  and  the  options  avail¬ 
able  are  described  in  the  mimeographed  pamphlet  "World  Map  Plot.  "  The 
mapping  function  remains  substantially  unchanged,  and  will  not  be  discussed 
here.  A  data  file  containing  the  latitude  and  longitude  of  the  necessary 
points  was  provided,  and  presently  exists  as  a  permanent  file 
(WMAP  DATAX2401,  CY=2,  ID=WINTERS).  This  must  be  attached  as  TAPE4, 
and  is  read  with  an  unformatted  read  in  MAPLOT. 
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The  main  program,  MAP,  has  six  sections.  These  are  an  initialization 
section,  four  optional  sections  which  generate  points  for  lines  to  be  drawn 
on  the  map,  and  a  final  part  which  calls  the  mapping  subroutines  to  do  the 
actual  drawing.  The  four  optional  sections  ate  controlled  by  the  integer 
parameters  NOGC,  DIP,  AUR,  and  NOL.l4,  which  are  read  at  statement  250 
(along  with  parameter  CRT).  A  value  of  zero  for  any  of  these,  causes  that 
section  to  be  skipped.  A  negative  value  permits — and  a  positive  value 
suppresses  —  the  printing  of  the  points  which  are  generated.  A  considerable 
amount  of  information  is  provided  in  the  comment  section  of  the  program 
listing.  In  particular,  there  is  a  description  of  the  use  of  each  input  para¬ 
meter.  In  some  cases,  the  result  of  selecting  each  of  the  possible  acceptable 
values  for  the  parameter  is  also  given  in  detail.  These  comments  are  found 
just  before  or  just  after  the  corresponding  read  statements.  The  user  should 
refer  to  them  for  detailed  information. 

2.  2  Great  Circles 

NOGC  is  the  number  of  great  circles  to  be  calculated  and  drawn.  The 
user  supplies  the  latitude  and  longitude  of  the  end  points  of  the  great  circle, 
and  the  program  produces  a  locus  of  points  along  it.  Optionally,  the  program 
can  be  told  to  use  the  latitude  and  longitude  coordinates  of  any  of  eight  pre¬ 
assigned  locations  as  one  end  point,  in  which  case  the  user  supplies  only  the 
other  end  point.  (The  eight  fixed  locations  are  the  VL.F  broadcast  stations 
of  the  worldwide  Omega  navigation  system.)  One,  several,  or  all  eight 
stations  can  be  automatically  used  as  one  end  point,  with  a  single  user- supplied 
point  for  the  other. 

This  section  requires  some  input  data.  The  necessary  parameters  are 
read  at  statement  100,  and  the  end  points  are  read  in  between  statements  120 
and  160.  A  detailed  explanation  is  given  in  the  program  listing. 
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On  a  spherical  earth,  a  great-circle  between  points  A  (0^,  and  B 
has  an  arc  length,  in  radians,  of 


♦  =  arccos  ^  sin  0  sin  0^ 


cos  0^  cos  0^ 


cos 


1) 


where  is  the  difference  in  the  longitude  coordinates  and  0  represents 
latitude.  Then  the  locus  of  points  (0,  $)  on  the  great  circle  is  given  by 


=  arc  sin  ^tan  a  tan  0^  -  arc  sin  ^tan  a  tan  0^  2) 

where  a  gives  the  arc  length  between  the  north  pole  and  that  point  on  the 
great  circle  nearest  the  pole: 


a 


arcsin 


cos  0  cos  0,  sin 
o  1 


sin  t 


3) 


Subroutine  DEC  is  used  to  calculate  the  arc-length  i)i.  Then,  using 
subroutine  POINTS  finds  the  locus  of  points,  which  is  the  desired  output. 
Subroutine  CONV  is  used,  at  the  beginning,  to  convert  the  end-point 
coordinates,  which  are  read  as  pairs  of  seven-digit  floating  point  numbers 
DDDMMSS  representing  degrees,  minutes,  and  seconds,  to  radians  and 
decimal  degrees. 


2.  3  Dip  Equator 

Any  nonzero  value  for  DIP  indicates  that  the  magnetic  dip  equator 
should  be  drawn.  The  latitude  points  are  stored  in  block  data  DIPEQ  and 
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in  CONtMON/C/  and  cor  respond  to  longitude  intervals  of  5°  starting  at 
-180°.  No  data  cards  are  required  for  this  section. 


2.  4  Auroral  Oval 

Any  nonzero  value  of  AUR  indicates  that  the  northern  auroral  oval 
should  be  drawn.  All  comments  pertaining  to  the  dip  equator  are  valid 
for  this  section  as  well. 


2.  5  Lines  to  be  Read  In 

NOLN  is  the  total  number  of  points  to  be  read  in,  if  other  lines  or 
contours  are  to  be  drawn  on  the  map.  The  coordinates  can  be  read  in 
-one  pair  per  card-  as  seven-digit  numbers  DDDMMSS  representing 
degrees,  minutes  and  seconds,  or  as  decimal  degrees.  In  the  former 
case,  subroutine  CONY  converts  them  to  decimal  degrees.  Details  are 
in  the  program  listing. 

If  more  than  one  distinct  line  is  to  be  read  in,  the  last  card  of  each 
line  must  contain  an  impossible  coordinate,  such  as  latitude  =  100.  This 
"break  point"  tells  the  plotter  to  lift  the  pen  when  going  from  one  line  to 
the  next.  These  points  are  counted  as  part  of  NOLN.  If  one  wishes  to 
avoid  counting  points  to  find  NOLN,  one  can  terminate  this  section  by 
setting  NOLN  equal  to  a  very  large  number  and  doing  one  of  two  things: 

a)  read  in  two  consecutive  cards  with  impossible  latitude  coordinates,  or 

b)  read  in  a  card  with  'X'  punched  in  column  1. 


2.  6  Storing  the  Points 

The  points  generated  in  the  optional  sections  are  all  found  before  any 
plotting  begins.  Initially,  they  are  stored  in  a  250x2  array  STORE.  Brea 
points  are  inserted  automatically  after  each  line  generated.  When  the 
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array  is  full,  the  numbers  are  written  to  a  scratch  file,  TAPE8,  along 
with  the  number,  NUM,  of  pairs.  The  array  then  is  filled  again,  and 
so  on.  At  the  end,  however  many  points  remain  in  STORE  are  written 
to  TAPE8.  After  the  map  is  complete,  TAPE8  is  read,  250  pairs  at  a 
time  and  the  lines  are  drawn  by  a  procedure  which  is  identical  to  that  for 
the  map  itself. 

2.  7  Miscellaneous  Comments 

This  program  hai  been  run  on  the  AFGL  CDC-6600  computer  using 
the  operating  system  SCOPE  3.4.4  and  the  FTN  compiler.  The  memory 
requirement  is  about  62000  octal  words,  slightly  above  the  system  default. 
The  compilation  time  is  about  7  seconds  and  a  single  map  requires  between 
10  and  15  seconds  beyond  this.  ‘Virtually  all  of  this  time  is  used  for  the 
map  itself;  very  little  is  used  for  the  great-circle  calculations. 

The  program  will  do  several  maps  consecutively.  A  full  set  of  data 
cards  must  be  read  in  for  each  map.  After  each  one,  the  program  returns 
to  statement  1  and  begins  reading  again.  A  blank  card  or  end-of-file  at 
that  point  terminates  program  execution. 

Parameter  CRT,  an  integer  read  at  statement  250,  should  be  zero  for 
PEN  plots  and  nonzero  for  CRT  plots. 

Whenever  (latitude,  longitude)  pairs  are  read  in,  a  field  of  20  spaces  of 
alphanumeric  information  beginning  in  column  22  is  also  read.  This  purely 
optional  feature  allows  one  to  identify  any  location  by  name,  and  have  the 
name  printed  on  the  output  listing  (not  the  map). 


3.  Examples 


3.  1  Use  of  the  Program 

The  following  deck  can  be  used  to  produce  a  map  with  the  on-line 
plotter 


WIN53,  CM65000.  3435  WIN  TER  STEINER 

REQUEST  (PLOT,*Q) 

ATTACH  (XX, .  file  containing  source 

FTN  (I=XX,  SL,  R=3)  program 

ATTACH  (PEN,  ONLINEPEN) 

ATTACH  (TAP E4,  WMAPDATAX2401 ,  ID=  WINTERS,  CY-2)  file  containing  map 
LIBRARY  (PEN)  data 

LGO. 

DISPOSE  (PLOT,  *OL) 

7/8/9 

data 

6/7/8/9 

The  following  list  contains  the  input  parameters  and  the  format  in  which 
they  are  to  be  used.  The  data  listed  are  sufficient  to  produce  a  view  of  the 
earth  from  above  the  North  Atlantic.  The  northern  auroral  oval  is  drawn,  as 
is  a  great  circle  from  Boston  to  Vienna,  and  another  line  (read  in)  connecting 
these  cities.  Great  circles  from  Bermuda  to  omega  stations  B,D,  F,  and  G 
are  drawn,  even  though  F  will  be  off  the  map.  The  map  will  be  8"  X  8”, 
drawn  with  pen  and  ink  with  grid  spacings  every  10°.  The  coordinates  of 
the  great  circles,  and  those  of  the  line  read  in,  will  be  printed.  Those  of 
the  auroral  oval  will  be  suppressed. 
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The  first  three  cards  are  the  map  parameters  (see  mimeographed 
pamphlet).  The  fourth  gives  the  program  options  (see  Section  II).  The 
next  three  give  the  parameters  and  end  points  for  the  first  great  circle, 
and  the  following  two  are  for  the  remaining  four  great  circles  (see  program 
listing).  The  last  set  of  cards  contains  the  coordinates,  in  decimal  degrees, 
of  a  single  desired  line. 


Orthographic  projection,  described  in  section  in.  1 


Miller  Modified  Mercator  projection  with  magnetic  equator, 
auroral  oval,  great  circles  and  one  line  read  in. 


Orthographic  projection,  view  from  above  North  Pole 
with  auroral  oval  and  great  circles. 
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VI.  PROGRAMS  FOR  PROCESSING  POLAR  VLF  PROPAGATION  DATA 


1 .  Introduction 

The  Propagation  Branch  of  the  Electromagnetic  Sciences  Division, 
R.  A.D.  C.  ,  is  engaged  in  ongoing  research  involving  continuous 
measurements  of  many  ionospheric  and  atmospheric  parameters. 

Analog  data  from  their  various  instruments  are  fed  to  a  data  acquisition 
system  which  digitizes  them  and  writes  them  onto  tape  for  storage  and 
subsequent  retrieval.  With  the  introduction  of  a  new  acquisition  system 
it  became  necessary  to  revise  the  data-handling  procedures.  Two 
programs  were  requested  and  developed  to  satisfy  the  requirements. 

The  objective  of  the  first  program  is  to  read  the  BCD- coded  field 
tapes  as  quickly  as  possible,  check  for  errors,  and  write  all  the  useful 
raw  data,  packed  as  efficiently  as  possible,  onto  a  backup  tape  in  the 
original  sequence.  The  objective  of  the  second  is  to  read  the  backup 
tape,  convert  the  raw  data  to  physical  units,  sort  it  efficiently,  and 
plot  it  channel  by  channel.  The  methods  by  which  this  is  done,  and  the 
subsidiary  functions  of  the  two  programs,  are  outlined  below. 

The  field  tapes  contain  information  from  many  different  sensors 
("channels").  A  single  scan  of  all  the  channels  is  identified  by  an 
"experiment  time."  The  field  tapes  consist  of  a  sequence  of  such 
scans,  each  of  which  includes  some  recognition  characters,  the  experi¬ 
ment  time,  some  "fixed  data”  which  remains  the  same  from  scan  to 
scan,  and  the  raw  experiment  data.  To  provide  some  organization  to 
this  stream  of  data,  a  certain  number  of  scans  are  grouped  into  single 
records,  which  can  be  up  to  4096  characters  in  length. 

In  the  present  arrangement,  40  channels  are  scanned  and  15  scans 
constitute  a  record.  The  scans  are  265  characters  in  length;  the 
records,  3975.  It  is  possible  to  prepare  field  tapes  in  other  formats, 
and  the  programs  accommodate  them  with  minimal  changes. 
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There  is  no  requirement  that  the  scans  be  equally  spaced  in  time, 
but  the  normal  circumstance  is  one  scan  every  30  seconds,  or  one 
record  every  7-1/2  minutes.  At  this  rate,  a  continuous  run  of  60  days' 
duration  would  produce  6,912,000  separate  data  points,  or  172,800 
points  per  channel,  in  11520  records. 

The  need  to  examine  such  large  amounts  of  information  is  the 
reason  for  the  plots.  Two  plot  options  are  offered:  1)  an  average 
over  a  certain  time  interval,  known  as  an  "averaging  period,"  and  2) 
the  scan-by-scan  results.  Even  in  the  latter  case,  time  is  organized 
into  averaging  periods  for  convenience. 


2.  Program  BACKUP 
A)  Reading  the  Field  Tape 

Complete  records  of  good  field  data  consist  of  15  identical  scans, 
in  the  present  arrangement.  The  BCD  code  for  each  of  the  265  charac¬ 
ters  in  a  scan  occupies  2  octal  digits,  or  6  bits.  When  a  record  is  read 
into  central  memory,  10  characters  fit  into  each  60  bit  word,  but  there 
is  no  particular  correlation  between  words  and  meaningful  groups  of 
characters  (such  as  single  data-values).  Appendix  B  illustrates  this 
point. 

The  assumed  format  for  each  scan  is 

"CLDDD<HH<  MM<  SSCL  TTTY  YJ±NNNNN±NNNNN . iNNNNN 

",  C,  L,  <,  and  (space)  represent  the  BCD  codes  for  those  characters. 
DDD  represents  the  BCD  code  for  three  digits  giving  the  day  of  the  year, 
while  HH,  MM,  and  SS  stand  for  the  hour,  minute,  and  second  of  the  scan 
in  question.  TTT,  YY,  and  J  stand  for  the  codes  for  digits  representing 


58 


1 


the  time  interval  (in  seconds)  between  consecutive  scans,  the  year 
(e.g.  78)  and  the  system  code  (1  or  2).  (Together  these  three  con¬ 
stitute  the  "fixed  data.") 

The  first  field  iNNNNN  stands  for  a  positive  or  negative  5 -digit 
integer  giving  the  data  in  the  first  channel,  and  there  are  40  such  fields 
per  scan. 

The  strategy  of  BACKUP  is  to  buffer  a  complete  record  into  ari  ay  A, 
and  interpret  it  with  one  DECODE  statement.  The  DECODE  operation 
translates  the  code  to  numbers  and  stores  them  in  15-element  arrays. 

For  example,  DDD  is  decoded  into  an  integer,  ddd,  and  placed  in  array 
T(NS,  1)  according  to  format  13.  (NS  identifies  the  scan,  1  to  15)  hh,  mm, 
and  ss  are  also  stored.  The  channel  data  are  put  into  STORE(NS,N), 
where  N  ranges  from  1  to  40.  The  special  symbol  "CL  is  stored  in 
SPS(NS)  and  the  other  special  characters  are  ignored. 

Any  records  with  parity  errors  are  ignored,  as  are  those  containing 
an  unexpected  number  of  words. 

B)  Checking  the  Records 

A  variety  of  abnormal  conditions  can  occur: 

1)  A  time  interval  other  than  normal  can  appear  between  con¬ 
secutive  scans  on  the  same  record,  or  on  different  records. 

2)  A  temporary  shutdown  after  some  scan  will  result  in  filling 
the  remaining  characters  of  the  record  with  the  BCD  code 
for  $. 

3)  An  extra  character  or  two  can  be  inserted  somewhere  in  an 
otherwise  normal  record. 

4)  The  fixed  data  on  the  tape  may  not  agree  with  the  expected 
fixed  data  values  which  are  read  in. 
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Note  that  if  2)  occurs,  1)  is  likely  to  occur  also.  Whenever  1)  occurs, 
the  program  makes  note  of  it  and  proceeds.  If  2)  occurs,  the  DECODE 
operation  has  failed  since  $  appears  where  a  number  shoidd  be.  However, 
the  program  is  prevented  from  terminating  or  printing  diagnostics  by 
previous  calls  to  ERRSET  and  SYSTEMC.  The  same  record  is  then  de¬ 
coded  scan  by  scan  until  the  first  error  is  encountered  again,  and  all  the 
good  scans  are  saved. 

Occurrence  of  3)  causes  the  DECODE  to  fail  also.  Scans  including 
and  following  the  one  with  the  first  extra  character  will  be  lost. 

Condition  4)  can  cause  rejection  of  the  record,  or  be  ignored,  at  the 
user's  option. 

After  the  various  checks  are  satisfied,  all  the  scans  of  useful  infor¬ 
mation  are  put  onto  the  backup  tape. 

C)  Packing  the  Backup  Tape 

For  each  field-tape  record  processed,  one  backup-tape  record  is 
written.  Array  B  is  used  to  store  the  information  as  it  is  processed. 

Each  backup  record  begins  with  a  single  60-bit  word  containing  the 
fixed  data  information  and  the  BCD  code  for  the  special  recognition 
symbol  which  requires  24  bits  (64646464  octal).  The  fixed  data 

are  represented  as  an  integer  tttyyj  and  use  the  other  36  bits.  The 
integer  30781  indicates  30-second  scan  intervals,  1978,  and  system  1. 

Each  of  the  15  scans  requires  15  words,  so  the  complete  record  is 
226  words  long.  The  first  of  these  15  words  contains  a  recognition 
symbol,  "  (6  bits)  followed  by  a  single  right- justified  integer  yydddhhmmss, 
such  as  78010023000  giving  the  experiment  time.  The  2nd  through  15th 
words  are  reserved  for  the  channel  data,  3  channels  to  a  word.  That  is, 
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the  first  20  bits  of  word  2  of  any  scan  contain  the  integer  corresponding 
to  NNNNN  from  channel  1;  the  second  20  bits,  from  channel  2;  and  so 
on.  The  last  40  bits  of  the  15th  word  of  each  scan  are  zeros,  and  a 
new  scan  begins  with  a  new  word. 

The  recognition  symbols  are  stored  as  BCD  codes  by  routine 
MXPUTX,  the  integers  as  simple  numbers  by  1SBYTX.  The  raw  channel 
data  range  from  -39999  to  +39999  (millivolts).  To  avoid  sign  problems, 
the  positive  data  are  biased  by  +40000  and  the  negative  data  are  stored 
as  their  absolute  values. 

A  buffer  out  operation  completes  the  processing  of  each  record. 

D)  Use  of  the  Program 
1.  Required  Input  Files 

a)  TAPE4:  The  field  tape,  which  should  have  an  L,  designation  on  the 
REQUEST  card. 

b)  INPUT  (TAPE5):  One  card-image  containing  the  following  information 


in  format  915,110.  Typical  values  are  given  in  parentheses. 

NFX  =  the  number  of  files  to  be  read  from  the  field  tape  (1) 

NWX  =  the  number  of  words  per  record  (398) 

NCH  =  the  number  of  channels  per  scan  (40) 

NSCAN  =  the  number  of  scans  per  record  (15) 

SCAN  =  the  expected  time  interval  between  scans  (seconds)  (30) 

YR  =  the  expected  year  in  which  the  experiment  began  (78) 

IC  =  the  expected  system  code  (1  or  2) 

NFLAG  =  0  or  1  (see  below)  (0) 

RMAX  =  total  number  of  records  to  be  read  (0) 

TMAX  =  last  experiment  time  to  be  processed  (dddhhmmss)  (0) 
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SCAN,  YR,  and  IC  are  to  be  compared  with  the  fixed  data  of  the  first 
scan  of  each  record.  If  NFLAG  =  0,  any  inconsistencies  will  be  noted 
but  otherwise  ignored.  If  NFLAG  =  1,  inconsistencies  result  in  rejection 
of  the  record.  If  RMAX  and  TMAX  are  zero,  the  program  reads  to  the 
end  of  the  last  file. 

2.  Memory  and  Time 

The  program  requires  about  41000  (octal)  words  of  core.  Execution 
time  is  approximately  .  1  sec  for  every  good  record  processed.  Frequent 
interruptions  increase  this  somewhat. 

3.  Miscellaneous 

The  time  on  the  backup  tape  is  stored  as  a  single  integer  yydddhhmms s 
The  program  recognizes  a  change  of  year,  going  for  example  from 
77365235959  to  78001000000  instead  of  to  77366000000.  It  also  recognizes 
leap  years. 

3.  Program  READ 
A)  Basic  Organization 

The  handling  and  plotting  of  data  are  controlled  by  a  table  of  directives, 
DRCTV,  which  is  read  in  and  which  contains  15  entries  for  each  of  the  40 
channels  which  could  be  used.  This  is  described  in  Section  HI-D  and 
APPENDIX  A. 

Of  all  the  data  found  on  the  backup  tape,  only  some  of  it  may  be  of 
interest.  A  "time  window"  giving  the  first  and  last  experiment  times 

desired  is  read  in,  and - except  for  the  time  words - only  scans  in 

this  range  are  read. 

The  sheer  volume  of  information  to  be  processed  makes  the  use  of  a 
direct-access  file  attractive.  Suppose  there  are  170,  000  scans  of  40 


channels,  to  be  read  in  the  original  time  sequence.  One  first  wishes  to 
plot  the  information  obtained  over  the  entire  time  range  for  just  one 
channel,  and  then  information  from  the  entire  time  range  for  the  next 
channel,  and  so  on.  Even  lumping  20  scans  to  an  averaging  period, 
core  storage  for  some  340,000  points  is  lacking.  Most  alternatives  to 
core  storage  involve  rewinding  and  rereading  one  or  more  sequential 
files  to  pick  out  every  40th  point  (or  even  every  40th  group  of  points) 
and  are  cumbersome  and  costly  in  time. 

The  WRIT  MS  and  READMS  mass  storage  routines  allow  one  to 
write  many  records  to  a  scratch  file  and  recover  them  in  any  desired 
order.  The  strategy  of  program  READ  is  therefore  to  read  the  tape 
until  a  large,  but  not  unmanageably  large,  number  of  points  has  been 
obtained.  These  are  sorted  according  to  channel,  and  when  the  array 
containing  them  (chosen  to  be  10240  points)  is  full,  separate  records 
are  written  containing  information  from  each  channel,  using  WRITMS 
on  the  random- access  file  TAPE8.  The  array  is  then  refilled  and  the 
process  repeats.  If  all  40  channels  are  used,  the  data  from  channel  1 
appeax  on  records  1,  41,  81,-...  which  can  be  recovered  in  order 
without  reference  to  records  2-40,  42-80,  etc.  Data  from  channel  2  is 
on  records  2,  42,  82,  ...  and  so  on. 

Note  that,  while  data  from  all  40  channels  is  on  the  backup  tape,  one 
can  disregard  as  many  of  the  40  as  one  wishes  after  reading  them.  Note 
also  that  the  "number  of  points  obtained"  depends  on  the  plot  option.  If 
averaged  data  are  plotted,  there  is  only  one  point  per  channel  per  aver¬ 
aging  period,  regardless  of  the  number  of  points  used  to  get  the  average. 
If  scan-by-scan  results  are  plotted,  all  data  points  are  "obtained." 

The  program  logic  is  best  understood  as  an  organization  of  the 
sequence  of  experiment  times  found  in  the  time  window  into  different 


blocks  of  time.  A  single  scan  occurs  at  an  instant  in  time.  Each 
averaging  period  contains  whatever  number  of  scans  that  fall  within  its 
boundaries.  (Averaging  periods  start  at  the  beginning  of  the  time  window 
and  run  consecutively  to  the  end;  if  there  are  time  gaps  on  the  tape,  some 
periods  may  have  fewer  scans  than  normal,  or  maybe  no  scans.  )  Each 
dump  contains  data  from  as  many  averaging  periods  as  can  be  stored  in 
10240  words.  As  such,  the  time  interval  represented  by  a  single  dump 
depends  on  the  number  of  channels  being  processed,  and  the  plot  option. 
(A  dump  is  actually  the  process  of  writing  records  to  TAPE 8.  ) 

No  more  than  1360  records  can  be  dumped  in  all,  so  if  there  is  more 
information  than  these  can  hold,  a  plot  must  be  performed  before  the 
backup  tape  can  be  read  any  further.  A  plot  period  is  whatever  time 
interval  it  takes  to  fill  1360  records.  In  normal  circumstances  it  is  not 
necessary  to  go  to  a  second  plot  period. 

For  sake  of  illustration,  suppose  that  only  32  channels  are  to  be 
plotted.  Then  each  dump  can  contain  320  averaging  periods,  since 
320x32  =  10240.  If  the  averaging  periods  are  600  seconds  in  length,  each 
dump  will  represent  192,000  seconds,  or  2.22  days.  At  32  records  per 
dump,  42  dumps  can  be  made  before  the  allowable  number  of  records  is 
exceeded.  A  plot  period  is  therefore  93.  33  days. 

Suppose  that  the  second  plot  option  is  desired:  every  scan  is  to  be 
plotted.  The  averaging  period  is  then  arbitrary  as  far  as  the  program 
function  is  concerned,  and  should  be  made  as  large  as  possible  for 
efficiency's  sake  (6 1  scans  per  period).  Then  in  the  example  above,  each 
dump  contains  320  scans,  which  represents  a  much  shorter  time  interval 
but  the  same  number  of  data  points  to  be  stored  and  plotted. 
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B)  Storage  and  Manipulation  of  Data 

The  main  subdivisions  of  the  program  are  READ,  subroutine 
REDUCE,  and  subroutine  PLMAT.  The  latter  has  entries  PLSAV, 
PLSND,  and  PLTT.  The  tape  is  read  in  READ,  the  data  are  reduced 
in  REDUCE,  and  data  are  written  to  mass  storage  in  PLSAV  and  PLSND. 
Later  on,  PLTT  does  the  plotting. 

Upon  being  read  from  the  tape,  data  are  moved  through  storage 
locations  as  follows: 

1.  READ:  The  tape  is  read  by  the  main  program,  READ.  Every  scan 
is  read  into  STOR,  a  40-element  array.  After  each  scam,  subroutine 
REDUCE  is  entered. 

2.  REDUCE:  So  long  as  the  current  averaging  period  is  still  being 
filled,  the  contents  of  STOR  are  simply  copied  into  columns  of  ARRAY 
(a  40  by  64  array)  starting  at  column  4.  Then  control  is  returned  to 
READ  for  a  new  scan.  An  averaging  period  can  have  up  to  61  scans. 

Once  a  scan  is  read  which  has  an  experiment  time  past  the  current 
averaging  period,  STOR  is  held  in  limbo  while  the  contents  of  ARRAY 
are  processed: 

i.  The  raw  data  in  each  row  of  ARRAY  are  searched  for  their 

minimum  and  maximum.  Their  average  value  is  computed,  and 
these  numbers  are  stored  in  the  first  three  elements  of  each  row 
of  ARRAY. 

ii.  The  raw  data  are  converted  to  physical  units  in  ARRAY. 

iii.  Most  of  ARRAY  is  copied  as  is  into  BRRAY,  which  is  also  40  by 
64.  Those  rows  of  ARRAY  corresponding  to  unused  channels 
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are  eliminated,  so  BRRAY  is  more  compact,  in  general. 

Columns  of  BRRAY  are  NREC  in  length.  (NREC  is  the  number 
of  records  per  dump,  since  there  is  one  record  per  useful 
channel  per  dump. )  This  transfer  eliminates  many  logical 
operations  later  on. 

iv.  Contents  of  BRRAY  are  printed,  if  desired. 

v.  PLSAV  is  called,  to  store  whatever  information  is  to  be  plotted, 
and  control  is  returned  to  REDUCE. 

At  this  point,  the  current  averaging  period  has  been  taken  care  of.  and 
indices  are  reset  for  the  next  one.  STOR  is  read  into  column  4  of  ARRAY, 
more  scans  are  read,  and  the  process  is  repeated. 

3.  PLSAV:  The  points  to  be  plotted  are  stored  in  CRRAY.  a  10240- element 
array,  prior  to  being  dumped  to  mass  storage.  CRRAY  is  implicitly 
partitioned  into  NREC  sections  (one  per  channel)  of  RECL  words  each. 

Each  section,  once  filled,  will  constitute  one  record  of  the  dump. 

For  plot  option  1,  there  is  one  point  per  channel  per  averaging  period. 
Section  1  of  CRRAY- --elements  1  through  RECL--- receives  points  from 
channel  1  for  averaging  periods  1  through  RECL;  section  2- -elements 
RECL+1  through  2*RECL— -receives  points  from  channel  2.  Section 
NREC,  the  last,  receives  points  from  the  last  channel  used.  For  aver¬ 
aging  period  NAP,  the  location  in  CRRAY  for  data  from  channel  N  is 
LOC  =  (N-1)*RECL  +  NAP.  That  is,  CRRAY(LOC)  contains  that  point. 

After  storing  one  point  per  channel  in  CRRAY,  control  returns  to  REDUCE 
unless  CRRAY  is  full. 


Once  RECL  averaging  periods  have  been  completed, 
of  CRRAY  is  full,  and  the  dump  is  made;  NREC  records 


the  useful  portion 
.  of  length  RECL, 
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numbered  consecutively  in  the  index  register  INDEX,  are  written  by 
calls  to  WRITMS.  For  dump  number  NDUMP,  the  index  of  the  record 
containing  data  from  channel  N  is  IX  =  (NDUMP-  1)*NREC  +  N. 


(Note  that,  prior  to  the  WRITMS  operation,  each  section  of  CRRAY 
is  copied  into  array  DUMP,  which  has  RECL  elements.  DUMP  shares 
memory  with  array  ARRAY  from  REDUCE,  which  contains  information 
that  is  superfluous  by  the  time  PLSAV  is  entered. ) 

After  a  dump,  control  returns  to  REDUCE,  and  the  new  averaging 
period  begins. 

For  plot  option  2,  there  are  many  points  per  channel  per  averaging 

period - say  20,  for  sake  of  illustration.  Then  section  1  of  CRRAY 

receives  20  points  from  channel  1  after  the  first  averaging  period,  20 
more  after  the  second,  and  so  on  with  the  same  true  of  all  sections.  If 
RECL  =  320  as  in  the  example  cited  earlier,  16  averaging  periods  would 
exactly  fill  CRRAY  and  a  dump  would  be  made.  If,  however,  there  were 
21  points  per  period,  then  after  the  16th  period  the  following  would  occur: 
the  first  5  scans  are  packed  into  each  section  of  CRRAY,  filling  them; 
the  dump  is  made;  and  then,  after  indices  are  reset,  CRRAY  is  refilled 
from  the  beginning  with  the  remaining  16  scans.  Then  control  returns 
to  REDUCE. 

To  summarize,  STOR  holds  data  from  one  scan;  ARRAY  and  BRRAY, 
from  one  averaging  period;  and  CRRAY,  from  one  dump.  The  records 
written  to  TAPE8  eventually  contain  all  the  information  to  be  plotted  in  one 
plot  period. 

Generally,  the  time  window  will  end  with  CRRAY  only  partly  filled. 
The  last  dump  will  then  produce  "short  records"  with  SRECL  words, 
rather  than  RECL  words.  This  is  accomplished  by  entering  PLSND 
instead  of  PLSAV. 
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C)  Plotting 

Once  the  full  time-window  has  been  read,  plotting  begins  in  entry 
PLTT.  An  array  must  be  established  for  the  time  axis,  which  is  the 

y-axis  on  paper.  For  plot  option  1- (Z)  the  array  contains - for  each 

point  in  the  first  record - the  decimal  number  of  days  (hours)  from 

the  beginning  of  the  time  window.  The  first  2560  elements  of  CRRAY 
(whose  original  purpose  has  been  served)  are  reserved  for  this. 

The  labels  for  the  time  axis- -which  must  be  drawn  NREC  times - 

are  of  the  form  DAY  161  DAY  162.  . .  or  JUN  1 0  JUN  11...  for  plot 
option  1;  or  900  1000  1100...  for  option  2.  These  are  stored  in 
CRRAY  also. 

The  time  axis  starts  at  the  beginning  of  the  first  day  (hour,  if 
option  2)  in  the  time  window.  Of  course  the  data  start  wherever  the 
time  window  starts. 

For  each  channel  the  complete  axes  are  drawn  and  labelled  first. 
Labels  other  than  for  the  time  axis  come  from  DRCTV.  Then  the  data 
records  are  retrieved  by  calls  to  READMS.  The  index  for  the  Ith 
record  corresponding  to  channel  N  is  LX  =  (I-1)*NREC  +  N.  Records 
are  read  into  array  DUMP,  which  constitutes  the  x- array  for  plotting. 

The  y-array  for  plotting  is  BRRAY,  which  must  be  constructed. 

For  the  first  record  of  a  given  channel,  it  is  identical  to  the  first  ele¬ 
ments  of  CRRAY;  for  succeeding  records  CRRAY  is  augmented  by 
however  many  days  (hours,  option  2)  correspond  to  the  number  of  full 
records  already  read  and  plotted.  In  other  words,  if  RECL  points 
correspond  to  DPR  days  (per  record)  then  the  y-array  for  the  Ith  record 
is  BRRAY (J)  =  CRRAY (J)  +  (I-1)*DPR,  for  J  =  l.RECL. 
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The  plot  routine  LINE  plots  DUMP  versus  BRRAY  for  the  first 
record  of  the  channel  in  question.  Then  the  second  record  is  read 
into  DUMP,  BRRAY  is  reconstructed,  and  the  next  RECL  pointB  are 
plotted.  This  process  continues  until  the  la6t  record  has  been 
retrieved  and  plotted.  Then  the  plot  origin  is  reset  and  the  next 
channel  is  processed. 

D)  Use  of  the  Program 
1.  Required  Input  Files 


a)  TAPE 3  contains: 

- A  table  of  coefficients  (COFS(N,  J),  J=  1 , 5)  for  reduction  of  VLF 

data  in  26  channels  (13  channels  of  phase  information,  13  channels 
of  amplitude  information;  N  =  1,  13).  13  card  images,  each  read 

in  FORMAT(5F10.  5) 


- A  table  of  integer  directives  (DRCTV(N,  J),  J  =  1 ,  15 )  for  plotting 

and  handling  data  for  each  channel  on  the  field  tape.  Appendix  A 
explains  the  table.  40  card-images  in  FORMAT  (213,  4A10,  313, 


315,  313). 

b)  TAPE4:  the  backup  tape  containing  the  experiment  data 

c)  INPUT  (TAPE5):  three  card-images  with  the  following  variables: 

- NFX,  NWX,  NSCAN,  NCHAN,  SCAN,  YR,IC  FORMAT(7I5) 

— -TBEG,  TEND,  AVG  FORMAT(2I10, 15) 

---PRINT,  PLLL,  ARAN,  FACTOR,  ifXP  SIZE  FORMAT(2I5,  2F5. 2,  215) 
Typical  values  are  given  in  parentheses.  Program  information: 

NFX  =  number  of  files  on  the  backup  tap\^  (1) 

NWX  =  number  of  words  per  record  on  the  backup  tape  (226) 
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NSCAN  =  number  of  scans  per  record  on  the  backup  tape  (15) 

NCHAN  =  number  of  channels  per  scan  on  the  backup  tape  (40) 

SCAN  =  number  of  seconds  between  consecutive  scans  (30) 

YR  =  year  in  which  experiment  began  (78) 

IC  =  system  code 

Time  window  and  averaging  period: 

TBEG  (dddhhmmss)  =  beginning  of  time  window 
TEND  (dddhhmmss)  =  end  of  the  time  window 


AVG  =  duration  of  a  single  averaging  period,  in  seconds  (600) 
Print  and  plot  options: 

PRINT  =  print  option:  0  means  no  data  are  to  be  printed;  1 

means  print  the  minimum,  maximum,  and  average 
value,  in  each  averaging  period;  2  means  print 
this  plus  the  scan-by-scan  results. 

PLLL  =  plot  option:  0  means  no  plot;  1  means  plot  the 
average  value  in  each  period;  2  means  plot 
scan-by-scan  results. 

ARAN  =  number  of  inches  per  day  (hour)  on  the  time  axis. 
DEFAULT  =  2" 

FACTOR  =  multiplicative  scaling  factor  for  plotting. 

DEFAULT  =  1.0 

IDT  =  code  for  labelling  time  axis:  0  means  "DAY  161 

DAY  162  1  means  "JUN  10  JUN  11  ..." 

DEFAULT  =  0 

SIZE  =  size  of  paper,  in  inches  (11"  or  30")  DEFAULT  =  30" 
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2. 


Other  Files  Used 


a)  OUTPUT  (TAPE6)  provides  running  commentary  on  anything 
unusual  on  the  backup  tape,  and  is  used  for  printed  data  output. 

b)  TAPE8:  random-access  scratch  file. 

c)  PLOT:  to  be  disposed  to  the  online  plotter. 

3.  Time  and  Memory 

The  program  requires  114,000  (octal)  words  of  core.  Processing 
1500  backup  tape  records  (about  8  days'  data)  takes  approximately  80 
seconds. 

4.  Miscellaneous 

Special  note  ought  to  be  made  of  the  plotting  procedures  if  PLLL  =  2. 
The  averaging  periods  are  just  bins  for  collecting  data  to  be  plotted.  No 
time  information  for  individual  scans  is  kept,  just  the  beginning  and  end 
of  the  period.  Therefore  a  problem  arises  if  there  are  time  gaps  on  the 
field  tapes.  All  data  within  an  averaging  period  are  plotted  from  the 
beginning  as  if  the  scans  were  equally  spaced  in  time,  and  consecutive. 

If  gaps  exist,  fewer  scans  occur  in  the  period.  The  averaging  period  is 
then  filled  at  the  end  with  minimum  values.  Therefore,  CAUTION:  if 
time  gaps  appear  on  the  field  tapes,  they  will,  in  general,  appear  at 
slightly  later  times  on  the  plot  than  they  actually  should.  Moreover,  if 
a  single  gap  runs  from  one  period  into  the  next  one,  it  will  appear  as 
two  gaps,  at  the  end  of  both  periods. 

If  this  becomes  a  problem,  the  averaging  period  can  be  shortened 
to  the  scan  time  (AVG  =  SCAN).  This  will  cause  the  program  exec  time 
to  increase  considerably. 
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If  PLLL  =  2,  the  number  of  scans  per  averaging  period  should  be 
constant.  Therefore,  AVG  Should  be  chosen  so  that  SCAN  divides  evenly 
into  it. 

The  preceding  comments  in  this  section  do  not  apply  if  PLLL  =  1. 

If  it  becomes  necessary  for  averaging  periods  to  accommodate  more 
than  61  scans,  the  dimensions  of  ARRAY  and  BRRAY  in  REDUCE  must  be 
increased.  The  present  40  by  64  size  allows  for  61  scans  and  the  minimum, 
maximum,  and  average  values  for  each  of  40  channels.  For  100  scans, 

40  by  103  would  be  needed.  Also,  LIMIT  in  subroutine  REDUCE  must  be 
increased,  and  the  dimensions  of  BRRAY  should  be  adjusted  in  PLMAT. 

If  more  than  40  channels  are  required,  it  is  sufficient  to  increase  the 
dimensions  of  any  arrays  which  are  presently  set  at  40,  and  augment  the 
data  reduction  section  of  REDUCE. 

Writing  and  reading  the  random-access  file  is  not  particularly  time- 
consuming.  A  quick  test  showed  that  500  records  of  20  words  cam  all  be 
created,  written  to  mass  storage,  and  retrieved  in  ’.7  seconds.  (For 
comparison,  if  these  records  are  also  output  to  the  printer,  the  WRITE 
operation  alone  requires  an  additional  1.5  seconds.)  Expanding  the 
record-length  to  500  words  only  triples  the  time  required  for  retrieval. 
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APPENDIX  A:  The  Table  of  Directives 

The  table  of  directives,  DRCTV  (see  section  Ul-D-l-a)  is  a  40  by  15 
array.  The  15  elements  of  each  row  have  the  following  meanings. 

1  (13)  Channel  index  (1-40).  Channels  are  numbered  0-39,  but 

indexed  1-40  because  of  the  FORTRAN  subscript  require¬ 
ment.  No  program  function,  only  used  to  keep  the  table 
in  order. 

2  (13)  Sensor  number  (1-40)  associated  with  a  given  type  of  instru¬ 

ment,  which  could  be  plugged  into  any  channel.  Used  for 
data  reduction  in  REDUCE. 

3  (A10)  10  characters  used  to  label  fee  data-axis  of  fee  plot. 

4  (A10)  10  characters  giving  fee  units  of  the  data-axis  label. 

5,6  (2A10)  20  characters  identifying  fee  type  of  measurement  being 
plotted. 

7  (13)  Data  disposition  directive:  0  means  ignore  all  the  data 

from  this  channel;  1  means  process  it.  In  fee  channel 
containing  y- component  magnetometer  data,  2  means 
compute  fee  horizontal  component  and  store  it  there. 

8  (13)  VLF  data  reduction  directive:  refers  program  to  correct 

elements  of  COFS. 

9  (13)  Reserved  for  riometer  subtraction  option. 

10  (15)  Maximum  data  value  to  be  plotted. 

11  (15)  Minimum  data  value  to  be  plotted. 

12  (15)  Total  number  of  inches  on  the  data-axis  for  this  channel. 

13-15  Not  used. 
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ARCON 


APPENDIX  B:  Dumps  of  a  Typical  Record 

The  following  is  a  character  dump  of  the  first  530  characters  (53 
words  corresponding  to  2  full  scans)  of  the  first  record  on  a  field  tape. 
Spaces  appear  in  the  record  exactly  as  shown. 


"C  L010<01<54<07CL  030782-00005+  9516+  0004+  9517+  0004+  9510+  0005+  9524+  0005+ 
9512+  0002+  9519+  0004+  5298+  0003+  5305+  0003-  4337+  0810-10007-  0012+  2991+  0 
016+  3038+  0757+  6724+  3296+  0003+  0002+  0002+  7367+  0003+  0002+  3856+  64?2-  005 
0+  3424-  0205-  0019+  1338"CL010<01<54<37CL  030782-00005+  9516+  0004+  95*7 +  0004 
+  9510+  0005+  9524+  0005+  9512+  0002+  9519+  0003+  5298+  0003+  5305+  0003-  4337+ 
0811-10017-  0012+  2882+  0016+  3038+  0766+  6652+  3201+  0003+  0002+  0002+  73^6+  00 
03+  0002+  4011+  6596-  0051+  3424-  0200-  0035+  1342 


The  following  is  an  octal  dump  of  the  first  31  words  (2  full  scans) 
of  the  backup-tape  record  which  was  made  from  the  field  tape  whose 
character  dump  is  given  above.  60-bit  words  are  separated  by  spaces 
for  clarity. 


64646464000000074076  64000001 105160457257  000001206026601l6l04  03013320470420140546 
02342120602720116105  03013200470410140557  023421005417101 16103  0260?6204704l401036l 
02373240116134000014  02474020470500124036  02371520555020124440  02342060470410116102 
02710160470414116102  02532400553040000062  02515000001464000023  02413640000000000000 
64000001105160457315  000001206026601l6l04  03013320470420140546  02342120602720116105 
03013200470410140557  02342020541710116103  02607620470414010361  02373260116134000014 
02474040470500124036  02371740554360124301  0234206047041 011 6102  02710140470414116102 
02537260554020000063  02515000001440000043  02413740000000000000 
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The  first  word  contains  the  recognition  character  64646464  plus  the 

octal  number  74076,  which  equals  30782  to  the  base  10  (the  fixed  data). 

In  the  second  word,  1105l60457257o  =  78010015407, „  which  is  the  first 

o  10 

experiment  time.  The  third  word  should  be  broken  down  into  bits  and 
reconstructed  in  thirds: 

00000120602660116104  equals 

8 

000000000000000001010000110000010110110000001001110001000100 

t  t  2 

The  first  20  bits  converted  to  decimal  give  5;  the  second  20  give  49516; 
and  the  third  20  bits  give  40004.  This  translates  to  -5,  +9516,  and  +4, 
which  can  be  seen  from  the  character  dump  to  be  the  raw  data  values  of 
the  first  three  channels  of  the  first  scan. 


VII.  IONOSPHERIC  PULSE  REFLECTIONS 

1 .  Introduction 

Substantial  modifications  were  made  to  an  existing  program  which 
had  been  written  to  calculate  the  effects  of  different  model  ionospheres 
on  VLF  signals  reflected  from  them.  These  effects  are  determined 
from  the  amplitudes  and  phases  of  the  four  complex  reflection 
coefficients,  |  I R  |  |  *  ||^j/  1^||'  and^R^*  The  first  subscripts 
refer  to  the  polarization  of  the  incident  field  with  respect  to  the  plane 
of  incidence  (the  vertical  plane  containing  the  wave-normal),  and  the 
second  refer  to  the  polarization  of  the  reflected  field.  In  an  aniso¬ 
tropic  medium  (S  f  0),  |  |R^  and  ^R  |  |  are  generally  not  zero  because 
of  Faraday  rotation  of  the  plane  of  polarization. 

In  the  case  at  hand,  a  fullwave  treatment  was  used  to  get  the 
phases  of  |  |R  |  |  and  |  |R^,  given  the  ground  distance  from  transmitter 
to  receiver,  as  a  function  of  the  angle  of  incidence.  For  nonuniform 
ionospheres,  these  are  found  by  successive  integration  of  four 
simultaneous  first-order  differential  equations  (relating  components 
of  the  electric  and  magnetic  field)  over  regions  in  which  the  ordinary 
and  extraordinary  indices  of  refraction  vary.  The  desired  quantities 
were  the  stationary  phases  of  |  |^  |  j  and  |  |R^ —  that  is,  the  angles  of 
incidence  at  which  the  phases  reach  their  maximum  values.  This  is 
sufficient  to  approximate  the  "path"  of  the  reflected  wave  to  a 
satisfactory  extent.  The  amplitudes  of  the  reflection  coefficients  at 
stationary  phase  and  the  group  heights  were  also  to  be  determined. 

2.  Program  Description 

The  program  is  written  so  that  a  range  of  frequencies  is  considered. 
For  each  frequency,  the  reflection  coefficients  are  calculated  (in 
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subroutine  BLBOX)  as  a  function  of  the  angle  of  incidence,  cp,  and 
the  phases  are  referred  to  some  specific  point,  usually  the  location 
of  the  transmitter.  The  phases,  originally  determined  modulo  2tt, 
are  adjusted  (in  subroutine  STAFAZE)  by  the  addition  of  multiples 
of  2rr  to  remove  artificial  discontinuities.  Then  a  term  is  added  to 
refer  the  phase  to  the  receiver  rather  than  the  transmitter.  The 
resulting  phases,  regarded  as  a  function  of  cp,  show  a  maximum 
somewhere  between  cp=Cf  and  cp=90°  • 

To  find  the  stationary  phase  angle  for  each  reflection  coefficient, 
the  table  of  phases  is  searched  for  its  greatest  value.  An  interpolation 
procedure  is  used  to  generate  a  new  table  including  only  the  region  of 
the  maximum.  This  in  turn  is  searched  for  its  maximum,  which 
should  be  very  close  to  the  desired  maximum.  As  a  check,  the 
reflection  coefficients  are  redetermined  at  three  points  straddling 
the  peak  and  the  actual  maximum  is  found  by  interpolation  using  these 
three  points.  This  procedure  locates  the  value  of  cp  corresponding  to 
stationary  phase  within,  at  least,  .05°.  The  values  of  the  amplitude, 
phase,  and  cp  for  each  reflection  coefficient  are  stored,  and  the 
program  proceeds  with  the  next  frequency. 

Occasionally  when  the  steps  in  cp  (PSTEP)  are  too  large,  consecu¬ 
tive  phases  actually  differ  by  more  than  2tt.  In  such  a  case,  an 
ambiguity  can  arise  in  STAFAZE  resulting  in  the  addition  of  the 
wrong  multiple  of  2tt.  The  program  recognizes  this  condition  and 
repeats  its  calculations  with  smaller  steps.  This  is  more  likely  to 
happen  for  higher  frequencies  than  for  lower  ones. 

The  range  of  cp  to  be  searched  for  the  maximum  in  the  phase 
angles  is  read  in  at  the  beginning.  If  the  maximum  lies  outside  this 


region,  the  program  prints  a  comment  to  that  effect  and  proceeds 
with  the  next  frequency. 

If  one  needs  the  values  of  the  reflection  coefficients  only  at  a 
single  angle  of  incidence,  one  can  set  PSTEP  equal  to  zero.  In 
this  case,  the  stationary  phase  section  is  bypassed. 

If  desired,  the  group  heights  for  the  two  cases  are  calculated 
in  subroutine  GRPHT,  after  the  stationary  phases  have  been 
determined  for  all  frequencies.  If  a  single  cp  is  used,  so  that 
stationary  phases  are  not  found,  this  section  is  bypassed. 

Another  available  option  is  to  plot  the  magnitude  of  the  reflection 
coefficients  at  stationary  phase,  and/or  the  group  heights,  as  a 
function  of  frequency.  If  stationary  phases  are  not  found,  the 
reflection  coefficients  generated  for  the  given  cp  are  plotted. 

It  can  happen,  for  model  ionospheres  with  profiles  containing 
particularly  large  gradients,  that  the  phases  do  not  increase  and 
decrease  monotonically  on  either  side  of  the  absolute  maximum. 

In  fact,  it  is  possible  to  get  "ledges"  (inflection  points)  or  even 
subsidiary  maxima  in  the  table  of  points  representing  the  phases. 

In  such  cases,  the  peak-location  procedure  can  give  erroneous 
results  or,  more  likely,  end  up  mistaking  this  problem  for  a  2tt 
ambiguity  somewhere.  In  the  latter  case,  the  program  cuts  the 
step  size  and  tries  again.  It  is  not  clear  whether  the  local  maxima 
are  due  to  actual  physical  circumstances,  or  rather  to  numerical 
procedures  in  the  supplied  routines  in  BEBOX  which  are  not 
suffid<mtly  accurate  for  the  sorts  of  ionospheric  models  encountered. 
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3.  Use  of  the  Program 

The  required  input  data  are  in  the  original  format,  except  that 
the  last  card  has  been  added  to  allow  the  group  height  calculations 
and  the  plotting  to  be  optional.  This  card  reads  three  integer 
parameters,  NPLOT,  GRP,  and  MPLOT  in  format  315.  If  NPLOT 
is  not  zero,  the  stationary  phase  plot  is  drawn.  If  MPLOT  is  not 
zero,  the  group  height  plot  is  done.  If  GRP  is  not  zero,  the  group 
heights  are  calculated.  If  any  of  these  parameters  is  zero,  the 
corresponding  section  is  skipped.  If  GRP=0,  MPLOT  is  automatically 
set  to  zero. 


A  list  of  input 

parameters  follows: 

Quantity 

F  ormat 

Explanation 

DG 

FI  0.  4 

transmitter-receiver  distance  (km) 

CMNT(1 ,  K),  K  =  1 , 6 

6A1 0 

general  comment 

CMNT (2,  K),  K  =  1 , 6 

6A10 

comment 

HNO,  DHN 

2F8.2 

starting  height,  height  interval 
for  electron  density  profile  (km) 

ENL(K),  K  =  1 ,  KN 

8F10.  5 

electron-density  profile 

CMNT (3,  K),  K  =  1 , 6 

6A1 0 

comment 

HNUO,  DHNU 

2F8.  2 

starting  height,  height  interval 
for  collision-frequency  profile  (km) 

ENUL(K),  K  =  1 ,  KNU 

8F1 0.  5 

collision  frequency  profile 

HREF 

FI 0.  2 

reference  height  (km)  (generally 
zero) 

HT(1 ),  NO(l ),  N1  (1 ) 

F10.  2,  215 

starting  height  for  integration  (km); 
number  of  integration  steps  per 
wavelength  (in  units  of  100);  step 
multiplier 
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HT  (2) 

F,  FH,  DIP,  AZ,  PHI 

FEND,  FSTEP,  PEND, 

NPLOT,  GRP,  MPLGT 


F10.2  stopping  height  for  Integration  (km) 

5F10.4  starting  frequency;  gyro  frequency; 

dip  angle;  azimuth;  starting  angle 
of  incidence 

PSTEP  4F10.4  last  frequency;  frequency  interval; 

last  angle;  angle  interval 

315  see  comments  above 
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VIII.  THE  IMPROVEMENT  FACTOR  FOR  THE  AR RESTED  SYNTHETIC 
APERTURE 

1 .  Introduction 

Program  ASAR  computes  the  improvement  factor  for  the  arrested 
synthetic  aperature  radar  (ASAR).  The  ASAR  processing  combines 
displaced  phase  center  antenna  (DPCA)  processing  with  doppler  filtering 
to  enhance  the  improvement  factor. 

2.  Problem  Background 

Suppose  an  aircraft  A  (velocity  v)  carrying  a  radar  antenna,  sees 
a  target  in  its  main  beam  at  an  angle  0^  (Fig.  1.  ).  At  the  same  time 
a  large  stationary  object  G  (ground)  appears  in  a  side  lobe  of  the 
antenna  pattern  at  angle  0.  A  portion  of  the  ground  G  has  an  apparent 
approach  speed  of  v  cos  0,  whereas,  a  stationary  target  in  the  mean 
beam  approaches  at  v  cos  0^.  It  is  possible  that  the  target  T  has  the 
same  apparent  speed  as  the  portion  of  ground  G.  Thus  conventional 
doppler  filtering  which  suppresses  the  main  beam  stationary  target 
cannot  suppress  side  lobe  ground  reflection.  One  method  of  suppress¬ 
ing  ground  clutter  is  ASAR  processing. 

3.  ASAR  Processing 

3.  1  General  Description 

In  ASAR  processing,  the  antenna  beam  used  for  target  detection 
is  assumed  to  be  formed  by  a  combination  of  two  antenna  arrays, 
array  A1  and  array  A2.  (See  Fig.  2).  Each  antenna  array  is  a 
linear  array  of  elements.  These  linear  arrays  are  aligned  along  the 
flight  path  (velocity  v).  Their  phase  centers  are  separated  by  a 
distance  d.  The  signal  received  by  each  antenna  array  is  processed 
by  identical  bandpass  filters.  The  outputs  of  both  filters  are  then 
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respectively  weighted  by  complex  constants  (processing  constant  A, 
and  processing  constant  B)  and  summed.  This  summed  coherent 
signal  is  then  detected  to  provide  an  appropriate  detection  statistic. 

The  signals  from  both  arrays  A1  and  A2  are  in  the  form  of 
pulses.  The  pulses  from  array  A2  are  emitted  at  a  slightly  later 
time  T  then  the  pulses  from  array  Al.  In  ASAR  processing,  the 
delay  time  t  is  chosen  to  be  the  ratio  between  the  separation 
distance  d  between  corresponding  array  elements  and  the  aircraft 
velocity  v  (t  =  d/v). 

Although  both  arrays  are  identical,  the  patterns  produced  by 
each  array  are  different.  This  difference  is  due  to  the  fact  that  the 
antennas  are  at  different  locations  in  the  airplane  and,  therefore, 
are  subjected  to  different  reflection  characteristics. 

By  appropriately  choosing  the  processing  constants  A  and  B, 
it  is  possible  to  partially  reduce  the  ground  clutter  (reflection  from 
ground)  in  comparison  to  the  signal  (reflection  from  target).  Indeed, 
if  the  patterns  from  both  arrays  were  the  same,  it  would  be  possible 
to  cancel  all  ground  clutter  by  choosing  the  processing  constant  B  to 
be  the  negative  of  processing  constant  A.  This  cancellation  of 
ground  clutter  is  due  to  the  fact  that  the  ground  is  stationary.  The 
target  signal,  however,  is  not  cancelled  because  it  is  moving. 

3.  2  The  Signal 

The  two  way  complex  voltage  (0)  received  by  the  i'th  array 
(i  =  1,2)  due  to  reflection  from  the  target  is  a  function  of  the  line  of 
sight  angle  0^  between  the  direction  of  flight  path  v  and  the  line  of 
sight  joining  the  aircraft  A  and  the  target  T. 
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The  net  complex  voltage  output  E  after  filtering  and  weighting  is 

a  weighted  coherent  sum  of  the  respective  complex  voltages 

(E  (9  ),  E  (0  ))  received  at  each  array,  provided  that  the  doppler 
1  o  c*  o 

frequency  is  within  the  frequency  band  (u)^  -  Aim,  u>c  +  Aiu)  of  the 
filter: 

r 

E  =  (  A  E,  (0  )  +  B  E_  (0  )eJU:oT,  uj  -  Au^.Suj  +  Au) 

|  1  °  20  c  d  c 
V  0  ,  otherwise 

The  doppler  frequency  U)^  is  the  doppler  frequency  associated  with 
the  radial  velocity  between  aircraft  and  target.  The  delay  time  T 
has  been  previously  defined. 

3.  3  Ground  Clutter 

In  addition  to  the  target  reflection,  there  is  ground  reflection 
(ground  clutter).  The  coherent  two  way  complex  voltage  E^  due  to 
an  increment  of  ground  at  a  line  of  sight  angle  6  is 

E  (0)  =  (A  E.  (6)  +  B  E  (0))  d0 
c  l  i. 

The  two  way  complex  voltage  received  by  different  increments  of 
ground  are  assumed  to  be  incoherent  with  respect  to  one  another. 
Therefore,  the  net  ground  clutter  power  C  is  the  sum  of  the  ground 
clutter  power  due  to  each  ground  increment. 
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If  the  clutter  power  C  is  normalized  to  the  power  from  a  single 
element,  we  obtain: 


C/C0 


tt 


6  +  A8 


0  -  A0 


c 


2 

|EC(9)|  d9 


The  clutter  angle  9^  is  the  line  of  sight  angle  that  corresponds  to  a 
doppler  frequency  Xc  equal  to  the  center  frequency  of  the  filter: 


4TTV 

A.  r  ,  ^ 

x  =  -  f„  sin  9 

c  c  R  c 


The  frequency  f^  is  the  carrier  frequency  of  the  radar.  The  angle 

interval  (6^  -  A8,  0c  +  A0)  corresponds  to  the  frequency  band 

(x  -  Ax,  x  +  Ax)  of  the  filter: 
c  c 


4tt  v 

A 

X  ±  Ax  =  -  f_.  sin  (0  ±  A9) 

c  c  R  c 


3 . 4  Noise 

The  noise  level  N.  from  each  antenna  is  the  sum  of  the  noise 

l 

levels  Nq  of  each  of  the  P  array  elements: 


The  antenna  noise  in  each  antenna  is  limited  by  the  doppler  band 
pass  filter  of  frequency  width  2  Au).  The  noise  output  N.  of  each 
bandpass  filter  is: 


N. 

l 


2Auj 


uu 


Prf 


N. 

l 


After  passing  through  the  band  pass  filter,  the  noise  hT  is  acted  on 
by  the  processing  constant  C.  (C  =  A,  C  =  B)  thus  producing  noise 

I  1  1  A 

N.: 

l 

Ni  =  I  ci  I2  Cj  =  A,  C2  =  B  . 


The  noise  produced  in  antenna  1  is  incoherent  with  respect  to 
the  noise  produced  in  antenna  2.  Therefore,  the  net  noise  output  N 
is  the  sum  of  the  noise  output  produced  by  the  first  antenna  and 
the  noise  output  produced  by  the  second  antenna: 

N  =  Nj  +  n!. 

It  is  convenient  to  express  the  above  net  noise  out  in  terms  of  the 
ratio  T]  of  the  noise  power  NQ  to  clutter  power  CQ  produced  by  one 
array  element: 


2^PTL  rA2  + 

uj  f 
prf  L 
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3.  5  Maximization  of  Improvement  Factor 

The  improvement  factor  I  is  the  ratio  of  the  signal  to  clutter 
plus  noise  ratio  S /  (C+N)  and  the  reference  signal  to  clutter  ratio 
Sq/Cq  due  to  one  array  element: 


I 


S/  (C+N) 

=  Vco 


The  signal  S,  clutter  and  noise  terms  N  are  positive  quadratic 
functions  in  the  processing  constants  A  and  B.  Therefore,  the 
improvement  factor  I  is  the  ratio  of  two  positive  quadratic  forms 
in  the  processing  constants  A  and  B: 


I 


C+  A  C 

+  __  ’ 
C  B  C 


It  is  desired  to  maximize  the  improvement  ratio  I  with  respect 
to  the  processing  constants  A  and  B.  If  the  positive  definite 
matrix  A  is  of  the  form  A  =  a  a  ,  and  the  matrix  B  is  positive 
definite  then  the  maximization  of  the  ratio  is  given  by: 


1 


a 
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4.  Program  ASAR 

Program  ASAR  computes: 

1.  The  improvement  factor  I  as  a  function  of  the  radial  target  velocity. 

2.  The  two  way  gain  functions  (0)  for  each  antenna  at  a  given  line 
of  sight  angle  0  vs.  the  doppler  frequency  cu^.  The  two  way  gain 
function  is  related  to  the  two  way  complex  voltage  as  follows: 

G.(0)  =  10  log  1  E.(0)  |2  i  =  1,2 

Plots  of  the  above  are  also  made. 

Program  ASAR  is  divided  into  the  following  parts: 

1.  A  subroutine  'PATTERN'  for  the  computation  of  the  two  way 
complex  voltage  pattern  due  to  a  linear  array  of  elements. 

2.  A  subroutine  ’IMPROVE1  for  the  computation  of  the  improvement 
factor. 
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IX.  COMPUTATION  OF  THE  SPECTRUM  OF  AN  ENSEMBLE  OF  SIGNALS 

1 .  Introduction 

Program  "SPEC1KM”  computes  the  spectrum  of  an  ensemble  of 
signals.  The  signals  are  constructed  from  the  data  obtained  from  a 
tape. 


2.  Reading  of  Tape 

The  tape  consists  of  a  large  number  of  files.  Each  file  consists 
of  several  records.  The  data  on  a  group  of  files  constitutes  the 
information  needed  to  construct  each  signal. 

3.  Unpacking  of  Data 

The  data  tape  is  a  9  track  tape.  One  track  is  for  parity  and 
8  tracks  are  for  storage  of  packed  words.  The  words  are  stored 
across  the  tape;  thus  each  packed  word  consists  of  8  bits. 

The  CDC  word  is  a  60  bit  word.  Fifteen  packed  words 
(characters)  are  read  into  a  pair  of  CDC  words.  Seven  packed  words 
plus  one  half  a  packed  word  are  contained  in  the  first  member  of  the 
pair  of  CDC  words.  The  remaining  half  of  the  last  packed  word  is 
contained  in  the  first  part  of  the  second  CDC  word. 

These  15  packed  words  that  are  stored  in  2  CDC  words  are  then 
unpacked  into  15  CDC  words.  Each  of  the  15  CDC  words  then  contains 
a  packed  word,  that  is  right  justified.  The  15  CDC  words  are  then 
divided  into  5  triplets.  Each  triplet  is  used  to  form  a  number.  The 
number  is  the  bit  string  formed  by  three  contiguous  packed  words. 

The  first  pair  of  CDC  words  at  the  beginning  of  a  file  contains  the 
identification  (year,  day,  test  no.  and  file  no.  )  for  the  file. 
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4.  Program  Spectrum 

Program  SPECTRM  is  divided  into  6  parts: 

1.  Reading  of  the  data  tape. 

2.  Unpacking  of  data  obtained  from  the  tape. 

3.  Construction  of  the  ensemble  of  signals. 

4.  Averaging  over  the  ensemble  of  signals  to  produce  one  averaged 
signal. 

5.  Computation  of  spectrum  of  the  averaged  signal  that  lies  within  a 
prescribed  time  window.  (The  spectrum  is  computed  using  the 
fast  Fourier  transform  method.  ) 

6.  Computation  of  amplitude  and  phase  of  spectral  components  of  the 
time  limited  averaged  signal. 

The  program  is  constructed  as  a  set  of  nested  loops.  The 
outermost  loop  governs  the  reading  in  of  files.  The  second  loop 
(going  inward)  governs  the  processing  of  a  record.  The  third  loop 
governs  the  unpacking  of  a  pair  of  words. 
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X.  MAXIMIZATION  OF  ANTENNA  GAIN 
1 .  Introduction 

Program  ANT  was  designed  for  the  maximization  of  antenna  gain 
in  curtain  arrays.  The  current  distribution  in  the  antenna  satisfies 
a  system  of  integral  equations.  This  set  of  integral  equations  is 
solved  using  a  matched  point  method  taking  the  classical  King  Mack 
Sandler  functions  as  a  basis. 

The  antenna  gain  is  a  function  of  the  current  distribution  and  is 
computed  in  the  standard  way.  The  current  distribution  is  a  function 
of  the  geometry  of  the  array.  Therefore,  the  antenna  gain  is  also  a 
function  of  antenna  geometry.  In  particular,  the  gain  depends  on  the 
element  lengths  and  inter-element  spacing. 

The  antenna  gain  may  be  maximized  with  respect  to  a  set  of 
parameters.  This  maximization  may  be  achieved  using  the  Hooke-Jeeves 
algorithm.  In  particular,  the  maximization  of  gain  may  be  achieved 
with  respect  to  element  lengths  and  inter- element  spacings. 


2.  Integral  Equation  for  the  Current  Distribution 
2.  1  Derivation 

The  current  distribution  in  an  antenna  may  be  derived  from 
Maxwell's  equation  relating  the  radiated  field-A,  to  the  current 
distribution  J: 


V  A(r,  t) 


1  5  A(r,t) 


=  -  4tt  J(r,  t) 


The  above  equation  may  be  Fourier  transformed  with  respect  to  time 
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to  give  the  following  equation: 


V  A  (r,  uj)  +  k  A  (r,  to) 


J  (r,  ou) 


(2) 


For  free  space,  the  above  differential  equation  may  be  converted 
to  Maxwell's  equation  in  integral  form: 


J 


A  =  I  G(r,  r')  J(r')  dr' 


The  free  space  Green's  function  is: 


G  (r,  r' )  =  e^/R 
R  =  I  r  -  r*  i 


(3) 


(4) 


If  the  antenna  consists  of  parallel  wires  oriented  along  the  z-axis 
(e.g.  curtain  array)  then  the  above  integral  (3)  may  be  decomposed 
into  integrals  along  each  wire: 


A  (r) 


G 


i  i  t 

(r ,  z. )  I.  (z. )  dz. 
J  J  J  J 


(5) 


The  index  j  indexes  the  elements  of  the  array. 


The  A-field  at  the  i'th  element,  is  given  by: 


i  t  i 

G„  {z,  z  )  I  (z  )  dz 


(6) 


The  Green's  function  G..  still  has  the  same  form  as  (4).  However, 
the  distance  R  is  the  distance  between  the  two  points  z  and  z 
located  on  element  i  and  j  respectively. 

The  form  of  the  A-field  may  also  be  deduced  from  a  combination 
of  symmetry  considerations  and  boundary  conditions.  Again,  if  the 
antenna  consist  of  elements  oriented  along  the  z-axis,  then  the  above 
Maxwell's  equation  (2)  in  differential  form  reduces  to: 


L 2  2 

+  r  |  A(z)  =  0 

Idz 


e  =  |k 


(7) 


The  above  equation  has  a  sinusoidal  solution  with  arbitrary  phase 
and  amplitude: 


A  (z)  =  a  cos  p  z  +  b  sin  {3  z 


(8) 


We  shall  assume  that  the  current  satisfies  the  following  symmetry 
condition: 


The  same  symmetry  must  be  inherited  by  the  A-field: 


A  (z)  =  A(-z) 


(10) 


This  symmetry  can  be  obtained  by  imposing  the  following 
restriction  on  the  above  solution: 


A  (z)  =  a  cos  p  z  +  b  sin  j3  |  z 


(ID 


Finally,  the  electromagnetic  field  must  satisfy  appropriate 
boundary  conditions.  Namely,  the  tangential  E-field  must  be 
continuous  at  the  wire  surfaces  of  the  antenna.  The  above 
boundary  condition  allows  the  A-field  at  the  i'th  wire  surface 
to  be  expressed  in  terms  of  the  gap  voltage  V  of  the  i'th  antenna: 


A.  (z) 


cos  P  z  +  -  Vk 


(12) 


The  intrinsic  impedance  of  the  surrounding  medium  is  T).  Equating 
the  two  expressions  (6,  12)  for  the  A-field  at  the  surface  of  the 
i'th  element,  we  obtain  the  system  of  integral  equations  for  the 
current  distribution  I.  (z). 
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— . 
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2.  2  Method  of  Moments 

The  method  of  moments  is  a  method  of  finding  an  approximate 
solution  to  a  system  of  integral  equations.  In  the  method  of  moments, 
the  current  is  expanded  in  terms  of  a  basis: 


I.  (z') 
J 


Jn(z) 


(13) 


The  system  of  integral  equations  (6,  12)  for  the  current  distribution 
then  becomes: 


E 


Q..  (z)  I.  =  A.  (z) 

ijn  jn  l 


I 


Q..  (z)  =  I  G  (z,  z  )  J  (z  )  dz 


n 


(14) 


We  may  choose  a  different  class  of  functions  (weighting 

functions)  that  also  span  the  same  space  as  the  basis  J  .  Taking 

n 

the  inner  product  on  both  sides  of  the  above  equation  (14)  with  the 
m'th  weighting  function,  we  obtain  the  matrix  equation: 


E 


<  W  ,  Q..  >  I.  =  <  W  , 

m  ijn  jn  m 


A.  (z)  > 


(15) 
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There  are  two  classes  of  weighting  functions  that  are  commonly 
used: 

1.  The  weighting  functions  that  are  the  basis  functions  (W  =  J  ) 

mm 

(Galerkin  method). 

2.  The  weighting  functions  that  are  Dirac  delta  functions  (matched 
point  or  collocation  method). 

In  the  matched  point  method  the  weighting  functions  are  chosen  to  be 
delta  functions: 


W 

m 


6  (z  -  zm)  . 


In  this  case  the  above  equation  (15)  becomes: 


(16) 


EQ..  (2  )  I.  =  A.  (z  ) 
ljn  m  jn  l  m 

j.n 


(17) 


We  see  that  in  the  matched  point  method,  both  sides  of  the 

equation  (17)  are  matched  at  a  sufficient  number  of  points 

so  as  to  be  able  to  determine  the  current  coefficients  I.  . 

Jn 


z 

m 
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A  good  set  of  basis  functions  are  the  King,  Mack  and  Sandler 
functions: 

Jj  =  sin  p  (h  -  |  z  | ) 

=  cos  p  z  -  cos  P  h 

=  cos  ^  p  zj  -  (cos  j  8  hj  (18) 

P  =  2tt/X 

h  =  height  of  dipole. 


For  an  antenna  in  free  space,  the  far  field  radiation  A  may  be 
determined  from  the  current  distribution  using  Maxwell's  equation 
in  integral  form  (5). 

Once  the  A-field  is  determined,  the  E-field  and  H-field  may  be 
simply  determined: 


E  =  jouA 

(19) 

H  =  -  V  x  A 
U) 


The  Poynting  vector  S  may  be  determined  from  the  E  and  H  field: 


S  =  E  x  H 


(20) 
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The  Poynting  vector  may  be  thought  of  as  the  amount  of  power  P 

A 

flowing  in  a  particular  direction  U. 


S  =  P  U  ,  P  =  |  S  |  .  (21) 

A 

The  power  flow  P  depends  on  the  direction  U.  The  direction  may 
be  specified  in  terms  of  spherical  angles  (9,  CP): 


U  =  cos  cp  sin  9  i  +  cos  $  cos  0  j  +  sin  §  k  •  (22) 


The  total  power  P  radiated  from  the  antenna  is  the  integral 
of  the  power  flow  P  (6,  cp)  over  all  directions: 


■I 


2tt  dcp  d0  sin  0  P  (0,  $ ) 


(23) 


The  antenna  pattern  G(0,  §)  (gain  pattern)  is  the  normalized  power 
flow: 


G(0,  cp)  =  P  (0,  *)/(Pt/4tt) 


(24) 


The  gain  G  (0,  cp)  in  a  particular  direction  (0,  $)  is  the  ratio  between 
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the  actual  power  flow  P  (0,  cp)  in  that  direction  and  the  power  P^,/4tt 
radiated  in  that  direction  by  an  equivalent  isotropic  antenna.  By 
an  equivalent  isotropic  antenna,  we  mean  an  isotropic  antenna  that 
radiates  the  same  total  power  as  the  antenna  in  question. 

A  parameter  that  characterizes  the  antenna  pattern  is  the  front 
to  back  ratio.  The  front  to  back  ratio  is  the  ratio  between  the 
maximum  gain  and  the  gain  in  the  opposite  direction. 

2.  3  Maximization  of  Gain 

The  antenna  gain  may  be  considered  to  be  a  function  of  a  certain 
set  of  variables.  In  particular,  the  gain  may  be  considered  to  be  a 
function  of  the  element  lengths  and  inter-element  spacings. 

There  are  several  methods  for  finding  the  maximum  of  a 
function  of  several  variables.  One  such  method  is  the  Hooke- Jeeves 
(HJ)  algorithm.  This  method  has  the  advantage  of  not  having  to  supply 
the  gradient  of  the  function.  The  HJ  method  is  a  systematic  searching 
method  for  finding  of  the  maximum.  The  search  begins  at  a  base 
point  bj.  A  local  search  (exploration)  is  established  about  the  base 
point  b^.  This  local  search  consists  of  changing  each  independent 
variable  x^  by  a  prescribed  step  size  Ek.  Each  variable  is  first 
increased  then  decreased  by  the  step  size.  If  a  better  point  is  found, 
it  is  adopted  as  the  new  base  point.  If  a  better  point  is  not  found, 
then  the  local  search  is  repeated  for  a  constricted  domain  (each 
prescribed  step  size  is  halved).  The  algorithm  is  terminated  when 
one  of  the  prescribed  step  sizes  falls  below  a  critical  value. 

Assuming  that  a  better  base  point  b^  is  established,  and  reasoning 
that  if  a  similar  exploration  were  conducted  about  base  point  b^  as 
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from  base  point  b^,  the  results  are  likely  to  be  the  same,  a  local 

search  about  the  second  base  point  b?  is  omitted.  Instead,  a  temporary 

base  point  t  is  established:  This  temporary  base  point  is  in  the 

direction  of  b,  to  b_,  but  of  twice  the  distance  from  b,  to  b  : 

12  12 


t 


2 


bl  +  2  (b2  -  bx) 


A  local  search  is  carried  out  about  the  temporary  base  with  the 
initial  step  sizes.  If  a  better  point  is  found  than  the  base  point  b^, 
this  point  is  established  as  the  third  base  point  b^.  We  then  proceed 
in  exactly  the  same  fashion  for  the  third  base  point  b^  as  for  the  first 
base  point  b^ . 

If  a  better  point  is  not  found,  then  a  retreat  back  to  the  second 
base  point  b^  is  made.  We  then  proceed  in  exactly  the  same  way 

for  base  point  b  as  for  the  first  base  point  b  . 

^  1 

2.  4  Program  ANT 

Program  ANT  is  a  program  that  is  based  on  the  matched  point 
method  for  solving  the  set  of  integral  equations  for  the  current 
distribution.  The  three  King  Mack  Sandler  functions  were  chosen 
as  a  basis.  Using  the  Hooke-Jeeves  algorithm,  the  gain  is  maximized 
with  respect  to  several  sets  of  independent  variables.  One  such  set  of 
independent  variables  is  the  set  of  element  lengths  and  inter-element 
spacings.  In  addition,  the  following  quantities  are  computed  at  the 
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maximum  point: 

1 .  The  driving  point  admittance  of  the  antenna. 

2.  The  relative  power  in  the  termination. 

3.  The  standing  wave  ratio. 

4.  The  absolute  directivity. 

5.  The  front  to  back  ratio. 


103 


10.  REFERENCES 


C.  J.  Drane,  "On  the  Three-Term  Theory  for  Analyzing  Thin 
Cylindrical  Dipole  Antennas,"  AFCRL-72-0744,  Dec.,  1972, 
Physical  Sciences  Research  Papers,  No.  523. 

Wilde  and  Beighller,  "Foundations  of  Optimization,  "  Prentice-Hall, 
1967. 

R.W.P.  King,  R.  B.  Mack  and  S.S.  Sandler,  "Arrays  of 
Cylindrical  Dipoles,"  University  Press,  1968. 


XI.  DETECTION  AND  TRACKING  OF  LOW  FLYING  TARGETS  BY 
GROUND  BASED  RADARS 

1 .  Introduction 

The  detection  and  tracking  of  low  flying  targets  by  ground  based 
radars  is  hampered  by  various  effects  such  as  terrain  screening, 
receiver  and  environmental  noise,  ground  clutter,  bird  clutter  and 
multipath.  A  computer  model  has  been  developed  to  simulate  the 
performance  of  a  ground  based  radar  so  as  to  include  the  effects  of 
terrain  screening,  bird  clutter,  ground  clutter,  Raleigh  distributed 
noise,  multipath  and  surface  roughness. 

The  analysis  of  the  tracking  simulation  proceeds  on  the  basis  of 
a  target  penetrating  the  radar  coverage  along  a  path  indicated  by  the 
dashed  line  on  Fig.  1.  The  path  is  identified  by  the  distance  of 
closest  approach  SD  and  the  target  is  located  along  the  path  by  its 
distance  from  the  point  of  nearest  approach.  As  the  target  penetrates 
from  the  right,  the  radar  range  equation  is  used  to  compute  the  signal 
to  noise  ratios  (S/N)  and  the  noise  to  clutter  ratios  (N/C).  These 
ratios  depend  upon  the  characteristics  of  the  radar,  the  electromagnetic 
scattering  cross  section  of  the  target  and  environmental  parameters 
(terrain  features,  weather,  noise,  electromagnetic  interference). 

2.  Analysis  of  Target  Detection  Procedure 

The  signal  to  noise  ratio  is  given  by  the  following  expression, 
which  may  be  derived  from  the  radar  range  equation 

(S/N)  =  PT  (X2)  GjJ  FUMF  (VR)  .  oT  (9)  F^L  (R^)/ ((4tt)2  •  Nq)  (1) 

where  P  is  maximum  power  of  compressed  pulse,  X  is  wavelength. 
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G  is  one  way  gain  of  transmit  receive  antenna,  FUMF  is  MTI  filter 
A 

function,  V  is  radial  velocity  of  target,  0  is  range  of  target  with 
R 

respect  to  y  axis, 0^,(0)  is  target  cross  section,  R^  is  range  of  target 

from  radar,  F.  is  function  for  specular  multipath,  N  is  total  noise 
ML  ° 

power  (Raleigh  distributed). 

The  noise  to  clutter  ratio  (N/C)  is  given  by 

N/C  =  (S/C)/ (S/N)  (2) 

where  signal  to  clutter  ratio  (S/C)  is  given  by 

<S/C>  =  4  <V  '  XmtI  <NH,0T(9)FUMF(V/[aMED-  VV  C'  T/2) 


Here  is  clutter  improvement  due  to  MTI  filter,  is  median 

cross  section  of  log-normal  clutter  per  unit  surface  area,  0^  is 
antenna  beamwidth  in  elevation  plane,  C  is  velocity  of  flight,  t  is 
compressed  pulse  width. 

The  probability  of  detection  P  ^  is  calculated  at  a  number  of 

points  on  the  target  trajectory.  The  calculations  are  performed  using 

3 

formulas  derived  by  Fante  . 

The  probability  of  false  alarm  (P^^)  is  the  probability  that  the 

received  signal  exceeds  some  threshold  level  Y  when  there  is  no 

3 

target  present.  Fante  has  shown  PFA  may  be  approximated  by  the 
expression: 
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where  erfc  (z)  is  complementary  error  function,  N/C  is  noise  to 

clutter  ratio,  NtT  is  number  integrated  per  scan,c  is  standard 
H  o 

deviation  of  clutter  cross  section  per  unit  surface  area. 

By  fixing  at  a  given  value  equation  (3)  is  solved  numerically 

for  the  threshold  Y  .  Equation  (3)  was  solved  using  Weg stein's 
4  ° 

iterative  method  . 

3 

Fante  has  also  shown  that  for  large  (S/N)  ratios  the  probability 
of  detection  P  of  a  target  with  log-normal  clutter  that  varies  from 
scan-to-scan  and  Raleigh  noise  that  varies  from  pulse-to-pulse  is 
given  by 


DT 


L  i  ) 

N  -1 

1  H 

Y 

o 

l  nh<s/n  )j 

|  exp 

1  +  Nh(S/N) 

(4) 


Computer  program  PDT  uses  equations  (1)  through  (4)  to  calculate 
P^rp  for  a  given  PpA  at  a  number  of  points  along  the  target's 
trajectory. 

3.  Simulation  of  Target  Tracking 

In  order  to  simulate  the  radar  tracking  of  aircraft,  it  is  necessary 
to  know:  (1)  the  probability  of  detection  as  a  function  of  radar 
parameters  and  environmental  conditions;  (2)  the  probability  that 
the  target  is  visible  due  to  terrain  screening  effects;  (3)  algorithms 
for  monitoring  target  track  together  with  appropriate  processing  to 
provide  accurate  estimates  of  target  position  and  speed  in  the  presence 
of  false  targets  (bird  clutter). 
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The  calculation  of  the  probability  of  detection  (Pqt^  *or  a  &*ven 
probability  of  false  alarm  Pp  includes  the  effects  of  target  fluctuations , 
noise,  clutter,  multipath  and  surface  roughness.  In  order  to  obtain 
realistic  simulations  of  target  tracking,  the  effects  of  terrain  screening 
must  also  be  included. 

In  the  process  of  radar  tracking  of  a  target,  two  sets  of  information 
are  of  value:  (1)  the  state  of  the  system  (target  position  and  velocity) 
and  (2)  the  future  state  of  the  system.  The  observations  of  the  system 
are  corrupted  by  noise.  The  optimum  estimate  of  a  set  of  observations 
may  be  obtained  by  passing  the  data  through  a  filter. 

The  simplest  filter  model  is  used  for  tracking,  i.  e.  the 
fading -memory  (exponential)  polynomial  filter.  The  effects  of 
association  and  correlation  are  neglected.  The  tracker  contains  a 
fixed  gain,  non-adaptive  a-(3  filter  and  initiates  track  on  two  hits  out 
of  three  looks.  Track  is  maintained  on  three  hits  out  of  five  looks. 

The  simplest  case  is  considered,  where  the  tracker  gives  no  report 
if  there  is  no  hit. 

The  computer  simulation  of  the  detection  and  tracking  of  aircraft 
proceeds  in  the  following  manner: 

1.  A  target  trajectory  is  considered  as  shown  in  Fig.  1.  For  given 

radar  parameters  and  given  environmental  conditions  the  probability 
of  detection  and  the  probability  of  target  visibility  P^^  are 

calculated  at  a  number  of  points  (R^,  Qp)  along  the  target  trajection. 

2.  Several  thousand  target  tracking  simulations  are  calculated.  In 
each  simulation  the  following  steps  are  performed: 
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(a)  a  uniformly  distributed  random  number  h ^  is  generated  for  each 
target  position. 


(b)  If  P 


pj  ’  Py  ^  s  a  hit  is  declared  at  position  (R  0  ); 

K  K 


otherwise,  the  look  is  declared  a  miss. 


(c)  Since  there  are  uncertainties  in  the  measurement  of  radar  range  R 

and  azimuth  0,  two  quantities  R  and  0  may  be  considered  as  Gaussian 

distributed  random  variables  with  mean  values  R  0  and  standard 

Is.  K 

deviations  AR,  The  standard  deviations  may  be  calculated  from 

the  formulas  given  by  Barton^. 


(d)  The  hit  or  miss  report  for  each  point  on  the  trajectory  together 

with  the  variables  RR,  0R  are  fed  into  the  a-p  tracking  filter.  The 

tracker  yields  at  each  point  on  the  trajectory  (R  0  )  a  report  with 

*  a  is.  i\  r 

estimates  (RR,  0R)  or,  no  report. 


T-i  A 

(e)  For  each  simulation,  the  absolute  errors  R  =  (R  -  R  ) 

E  K  K 

®  ~  are  calcu^ated.  Errors  in  speed  are  also  calculated. 


3.  These  errors  are  appropriately  summed  over  the  number  of 
simulations  and  divided  by  the  number  of  track  reports  to  find  their 
corresponding  mean  values.  If  there  is  no  track  report  in  a  given 
simulation,  the  errors  are  assumed  to  be  zero.  Also  the  variances 
in  the  corresponding  errors  are  calculated. 


4.  Results  and  Conclusions 

On  Fig.  3,  the  effects  of  multipath  are  suppressed.  The  curve 

labeled  no  clutter  corresponds  to  aMED  =  -54  dB  and  the  curve  labeled 

clutter  corresponds  to  =  -34  dB. 

MED 
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It  may  be  noted  how  the  high  clutter  case  seriously  reduced  the 

probability  of  detection,  even  though  the  probability  of  false  alarm 

-5 

is  relatively  high  (P  =  10  )  and  the  MTI  clutter  improvement 

factor  is  48  dB. 

There  are  conditions  when  the  radar  signal  reflected  from  the 
target  along  the  direct  path  may  interfere  with  energy  that  assumes 
from  the  ray  reflected  from  the  ground  (see  Fig.  4).  This  can 
produce  areas  of  varying  signal  strength  which  lead  to  regions  of 
varying  .  The  regions  of  low  P  can  cause  missed  target 
detections  and  loss  of  track.  The  detections  effects  of  multipath 
depend  upon  antenna  height,  target  height,  target  range,  dielectric 
properties  of  the  ground,  surface  roughness,  the  antenna  pattern  in 
the  elevation  plane  and  the  wavelength.  When  the  surface  from  which 
the  indirect  energy  is  reflected  becomes  rough,  the  specular  multi- 
path  interference  effect  is  decreased.  However,  diffuse  multipath, 
which  may  look  like  noise,  will  increase  as  the  reflecting  surface 
becomes  rougher. 

On  Fig.  5,  the  clutter  was  low.  The  curve  labeled  "rough  surface" 
corresponds  to  rms  height  of  surface  irregularities  ff  =  5m  and  the 

XT 

curve  labeled  "smooth"  corresponds  to  o  =  0. 

XT 

On  Fig.  6,  the  decrease  in  P  .p  due  to  the  combined  effect  of 

clutter  and  multipath  is  illustrated.  The  curve  labeled  "no  clutter" 

corresponds  to  n  =  -54  dB  and  multipath  effects  are  suppressed. 

The  curve  labeled  "clutter"  and  multipath  corresponds  to  =  34  dB. 

ME  D 

It  may  be  noted  that  multipath  improver  performance  at  some  ranges 
but  degrades  it  in  others.  The  effects  of  multipath  combined  with 
clutter  can  make  it  difficult  to  establish  and  maintain  target  track 
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and  discriminate  between  real  targets  and  bird  clutter. 

On  Fig.  7,  surface  is  smooth  and  =  34  dB.  The  difference 

MED 

in  ^-yp  when  the  probability  of  false  alarm  varies  from  10  to  10~ 
is  illustrated  on  Fig.  7.  The  primary  effect  of  clutter  is  to  require 
that  the  detection  threshold  be  raised  in  order  to  keep  P  within 
reasonable  bounds.  This  demonstrates  that  the  system  will  have  to 
operate  at  false  alarm  rates  higher  than  usual  in  order  to  maintain 
reasonable  values  of  ,  thus  requiring  particular  attention  to  be 

given  to  design  of  the  tracker. 
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